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ABSTRACT^X  v 

A  study  of  the  flight  mechanics  of  a  three  surfaced 
fighter  aircraft  was  accomplished  using  bifurcation  and 
catastrophe  theories.  The  emphasis  of  this  /••research  was 
to  find  or'J&evelop  nonlinear  systems  analysis  techniques 
capable  of  analysing  and  predicting  the  complex  behavior 
found  in  high  -  angle  —of  —  attack  flight.-.  A  primary 
selection  criteria  for  the  method  to  be)  used,  in  the 
research  was  that  the  method  be  more  gJrdoal  in  nature 
than  the  techniques  currently- — in  use7  The  technique 
selected  was  based  on  thc/concepts  found  in  bifurcation 
and  catastrophe  theories,  jfi he  primary  Method -of> research  ’oir  jla)  ) 
involved  <5&e^  calculation  of  equilibrium  surfaces ^where 
the  time  rate  of  change  of  each  of  the  state  variables  is 
set  to  zero  in  the  aircraft  equations  of  motion,  and  the 
states  that  will  satisfy  that  condition  are  found  for 
various  control  deflections.  The  equilibrium  surfaces 
were  calculated  numerically  through  use  of  continuation 

>An  equilibrium/  trajectory  analysis  of  a  fifth  order 
model  of  the  airpxaft  showed  significant  adverse  yaw 
problems  at  high /angles  of  atteefe  as  well  as  the  presence 
of  roll  coupling,  ‘phonemena.-^  Bifurcational  behavior  was 
found  in  the  same  fifth  order  model  and  verified  by  i 

research# ^previously  accomplished A  mechanism  for  <£ba»ZE 
transfer  of  kinetic  energy  during  the  catastrophic 

behavior  shown  is  hypothesised.  A  spin  equilibriuim 

surface  for  the  aircraft  was  determined  from  an  eighth 
order  model  the — ai reraft, and^  from  that  surface, 

possible  control  deflections  for  recovery  from  the  spin 
were  obtained,  Numerical  simulations  of  the  model  were 
accomplished,  out  few  useful  results  were  obtained  due  to 

the  unstable  n  iture  of  the  aircraft.  A  number  of  areas 
for  improvement  of  the  aircraft  design  and  methods  used 
as  well  as  poisible  areas  for  further  research  are 
identified. 
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CHAPTER  I 


INTRODUCTION 

A.  Historical  Perspective 

The  prediction  of  an  aircraft's  dynamics  and 
stability  behavior  has  been  a  topic  of  interest  since 
heavier  than  air  flight  was  first  conceived.  From  the 
very  beginning,  aircraft  designers  were  at  odds  in  their 
beliefs  of  how  aircraft  should  be  designed  in  terms  of 
stability.  For  example,  the  Wright  Brothers  favored  a 
slightly  unstable  design  to  increase  manueverability 
while  the  European  aircraft  design  community  opted  for 
much  more  stable  designs  [1].  The  Wright  Brothers 
insistence  on  maintaining  the  elevator  forward  of  the 
wing  in  their  early  designs  was  motivated  by  both  the 
desire  to  make  the  aircraft  a  servant  to  the  controlling 
desires  of  the  pilot  and  the  desire  to  provide  a  pitch  up 
moment  following  a  stall  to  prevent  the  dive  that 
typically  followed  stalls  in  aircraft  of  that  period. 
It  was  such  a  stall  initiated  dive  that  caused  the  death 
in  1896  of  Otto  Lilienthal,  the  most  prominent  of  the 
European  airplane  designers  at  the  time  [2], 

Being  bicycle  designers  and  riders,  the  Wright 
Brothers'  concept  of  controlling  an  aircraft  in  flight 
was  much  the  same  as  that  of  riding  a  bicycle. 
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The  bicycle  is  inherently  unstable  but  the  instability  is 
easily  overcome  by  the  rider  through  the  use  of  the 
handlebars.  Unfortunately,  as  aircraft  became  more 
sophisticated,  the  large  demands  placed  on  a  pilot 
attempting  to  control  an  aircraft  with  the  tail  forward 
of  the  wings  made  the  design  impractical  and  often 
dangerous.  The  tail  forward  design,  later  known  as  a 
canard  due  to  its  resemblance  to  a  duck,  was  to  be 
shunned  by  the  worldwide  aircraft  design  community  for 
many  years  to  come  [2].  It  was  during  this  period  that 
the  desire  for  static  stability  became  an  unquestioned 
design  feature  for  nearly  sixty  years.  As  will  be 
discussed  below,  the  Wright  Brothers'  understanding  of 
the  relationship  between  stability  and  manueverability 
was  far  ahead  of  their  time.  They  understood  that 
stability  and  manueverabilty  are  relative  terms,  and  that 
they  describe  positions  along  the  same  spectrum  with 
stability  being  at  one  end  of  the  spectrum  and 
manueverability  at  the  other.  Consequently,  you  can  not 
enhance  one  without  sacrificing  the  other.  The  most 
simple  definition  of  stability  when  applied  to  an 
aircraft  is:  a  characteristic  of  an  aircraft  in  flight 
that  causes  it  to  return  to  a  position  of  steady  flight 
after  being  disturbed  from  that  position. 
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Nearly  all  aircraft  possess  enough  stability  to  recover 
their  equilibrium  after  small  disturbances;  however,  none 
possess  the  stability  to  recover  its  equilibrium  under 
all  conditions.  The  need  to  know  the  conditions  under 
which  an  aircraft  will  maintain  or  lose  its  stability 
form  the  motivation  for  this  study. 

B.  Aircraft  Design  Considerations 
The  capabilities  of  a  combat  aircraft  are  often 
measured  in  terms  of  performance  factors  such  as:  wing 
loading,  maximum  thrust  to  weight  ratio,  engine  response 
time,  maximum  turn  rate  and  rate  of  turn  build  up, 
maximum  transient  and  sustained  load  factor,  minimum  turn 
radius,  maximum  and  minumum  specific  excess  power.  In 
addition  to  these  there  is  still  one  other  very  important 
measure  of  a  combat  aircraft's  worth.  This  performance 
measure  is  aircraft  agility.  Agility  is  the  ability  to 
change  the  direction  and  magnitude  of  the  aircraft 
velocity  vector  in  three  dimensional  space  and  to 
accomplish  this  more  rapidly  than  an  opponent  [3].  A 
measure  of  agility  could  be  determined,  in  a  sense,  by  a 
formula  or  set  of  formula  that  combines  each  of  the 
individual  factors  listed  above. 
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Aircraft  performance  in  terms  of  the  factors 
discussed  above  has  been  steadily  improved  with  each  new 
generation  of  combat  aircraft  so  that  today's  modern 
combat  aircraft  are  much  more  agile  than  their 
predecessors.  This  level  of  increased  agility  is  a 
direct  result  of  improvements  in  the  state  of  the  art  in 
each  of  the  following  areas:  stability  augumentation, 
structural  technology,  power  plant  capabilities  and 
aerodynamic  design.  A  considerable  increase  in  agility 
has  been  made  by  simply  returning  to  the  farsighted 
design  philosophy  of  the  Wright  Brothers.  For  example, 
like  the  Wright  Flyer,  modern  fighter  aircraft  now  use 
statically  unstable  designs  to  enhance  manueverability. 
This  concept  is  demonstrated  in  the  designs  of  the  F-16, 
F-17  and  F-18  fighter  aircraft.  These  aircraft  are 
statically  unstable  but  are  kept  under  control  by 
computerised  flight  control  systems  so  that  the 
manueverability  advantages  of  static  instability  can  be 
realised  without  overburdening  the  pilot.  Modern, 
lightweight  composite  materials  and  aircraft  structural 
design  improvements  have  made  it  possible  to  decrease 
aircraft  weight  while  at  the  same  time  increasing 
structural  strength.  This  produces  relatively  light 
weight  aircraft  that  have  load  factor  capabilities  of  up 
to  12  g* 8  thereby  giving  the  aircraft  distinct  increases 
in  agility. 


The  powerplants  of  modern  fighter  aircraft  are  now 
powerful  enough  to  take  advantage  of  both  the  computer 
aided  stability  systems  (cass)  and  the  increased  load 
factor  capabilities  of  our  current  fighter  aircraft. 
This  makes  it  possible  to  sustain  very  high  turn  rates 
and  to  sustain  flight  at  very  high  angles  of  attack. 
The  modern  powerplants  also  can  achieve  thrust  to  weight 
ratios  in  excess  of  one  thus  providing  substantial 
advantages  in  air  combat  situations.  In  the  area  of 
aerodynamic  design,  however,  progress  of  the  same 
magnitude  as  in  the  three  areas  discussed  above  has  not 
been  made.  Some  innovative  uses  of  the  control  surfaces 
along  with  increased  use  of  some  of  the  powerful  vortex 
interactions  around  the  aircraft  have  been  accomplished, 
but  for  the  most  part  little  has  been  done  to  enhance  our 
understanding  or  our  ability  to  -analyze  the  extremely 
complex  aerodynamic  behavior  of  an  aircraft  in  air 
combat  type  situations. 

A  possible  aircraft  to  fill  the  need  mentioned  above 
is  being  designed  and  tested  by  McDonnell  Douglas  and  the 
Air  Force  Wright  Aeronautical  Laboratory  at  Wright.  - 
Patterson  AFB,  Ohio.  The  aircraft  is  currently  known  as 
the  F-15  STOL  Demonstrator.  It  is,  basically,  an  F-lSB 
aircraft  with  a  close  coupled  canard  mounted  just  forward 
of  the  wing  on  the  engine  inlet  shelves. 
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Provisions  are  also  being  made  to  implement  thrust 
vectoring  engine  nozzles  (to  enhance  longitudnal  control) 
and  thrust  reversal  capabilities.  The  potential  benefits 
of  such  a  configuration  are  well  documented  [4,5].  A 
few  will  be  briefly  discussed  here  for  the  sake  of 
completeness.  The  close  coupled  canard  will  improve  the 
aircraft's  lift  in  the  high  angle  of  attack,  (aoa),  range 
and  increase  stalling  aoa.  This  is  due  primarily  to  the 
vortex  interaction  between  the  canard  and  the  wing. 
Regarding  the  canard  as  an  extra  lifting  surface  gives 
the  capability  of  direct  force  manuevering.  Direct  foroe 
manuevering  allows  the  manuevering  modes  to  be  decoupled 
from  each  other  so  that  turns  without  using  bank  and 
altitude  changes  without  changing  pitch  etc.  are 
possible.  Another  positive  effect  of  the  canard  is  that 
it  provides  the  ability  to  decrease  trim  drag  at  all 
speeds.  The  extra  lifting  surface  area  of  the  canard 
alone  offers  some  distinct  advantages  over  a  similar 
aircraft  without  a  canard  by  decreasing  overall 
wingloading.  The  extra  control  surface  also  increases  the 
possibilities  for  control  reconfiguration  enabling  a 
safe  return  to  base  following  battle  damage.  The  canard 
configuration  also  dramatically  increases  pitch  rates  and 
turn  rates  as  well  as  the  ability  to  obtain  those  high 
rates  quickly. 
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The  F-15  STOL  Demonstrator  project  is  still  in  the 
testing  and  development  stage.  Large  amounts  of  wind 
tunnel  testing  (static  and  forced  oscillation)  has  been 
accomplished.  A  reasonably  accurate  model  of  the  aircraft 
has  been  developed  through  use  of  the  data  obtained  from 
the  wind  tunnel  testing.  Full  scale  computerised 
simulations  as  well  as  'real  time*  pilot  in  the  loop 
simulator  testing  in  a  fully  motion  based  simulator  are 
also  currently  being  accomplished.  Actual  test  flights 
of  the  aircraft  are  scheduled  to  occur  in  March  1988. 


Originally,  aircraft  design  and  testing  involved 
more  trial  and  error  techniques  than  actual  engineering 
analysis.  Experimental  data  was  not  only  unreliable  but 
was  also  hard  to  acquire.  Therefore,  most  of  the 
aircraft  development  process  relied  on  a  'lets  try  it  out 
and  see  what  happens*  attitude.  As  aircraft  designs 
increased  in  complexity,  it  became  not  only  more  and  more 
expensive  to  design  aircraft  in  that  manner  but  also 
quite  dangerous.  It  was  apparent  that  more  work  had  to 
be  accomplished  on  the  ground  before  testing  the  design 
in  flight.  High  aoa  stability  and  control  studies  during 
the  forties  and  fifties  concentrated  mainly  on  obtaining 
a  qualitative  description  and  understanding  of  the  flight 
dynamics  observed. 
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A  lack  of  reliable  wind  tunnel  and  flight  testing 
techniques,  combined  with  little  or  no  computational 
capabilities  stood  in  the  way  of  any  major  progress  in 
the  area  of  high  aoa  flight  dynamics.  Conventional 
analysis  techniques  based  on  several  unrealistic 
assumptions  and  linearisation  were  used  to  analytically 
predict  aircraft  behavior  [6,7]. 

Experimental  techniques  and  facilities  as  well  as 
computational  methods  and  capabilities  have  undergone 
major  improvements  in  the  past  years.  Mow,  quite 
accurate  aircraft  models  obtained  through  modern  testing 
techniques  are  used  in  full  six  degrees  of  freedom  (dof ) 
simulations  with  more  accurate  results.  Additionally,  a 
large  number  of  advanced  mathmatical  techniques  have  been 
developed  to  deal  with  specific  nonlinear  phenomenon. 
These  techniques  when  combined  with  information  gained 
from  various  numerical  simulations  have  given  some 
insight  to  a  few  particular  high  aoa  flight  phenomena, 
but  there  is  a  conspicuous  lack  of  global  analysis 
techniques.  The  highly  nonlinear  and  very  complex 
behavior  of  an  airoraft  during  high  aoa  flight  coupled 
with  poorly  understood  aerodynamics  has  created  a 
tremendous  amount  of  resistance  to  progress  in  the  area 
of  developing  general  or  global  analysis  techniques. 
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The  development  of  a  global  analysis  technique  will 
be  of  great  benefit  to  understanding  and  possibly 
conquering  high  aoa  flight.  A  global  type  analysis  could 
offer  specific  guidance  as  to  when  use  of  a  particular 
advanced  mathmatical  technique  would  be  practical.  A 
more  global  picture  of  the  aircraft  dynamics  and  behavior 
will  enable  engineers  to  detect  areas  where  problems  are 
likely  to  occur  and  study  those  areas  more  closely.  The 
heavy  use  of  computer  simulations,  while  accurate  within 
the  area  investigated,  still  tends  to  give  a  * tunnel 
vision*  picture  of  the  aircraft  dynamics.  The 
possibility  of  missing  some  very  dangerous  or 
catastrophic  behavior  is  always  present  in  such  an 
analysis.  The  more  global  the  analysis,  the  less  the 
possibility  of  having  a  very  costly  surprise  during 
actual  flight  testing.  In  spite  of  the  great 
difficulties  discussed  above,  the  need  for  a  more  global 
approach  to  studying  high  aoa  flight  dynamics  is  stronger 
than  ever.  Modern  fighter  aircraft  currently  in  the 
inventory  are  able  to  operate  routinely  in  areas  that 
were  once  thought  of  as  outside  of  the  flight  envelope. 
As  the  ability  to  generate  extreme  pitch  rates  and  to 
sustain  very  high  angles  of  attack  continues  to  improve, 
we  must  also  improve  our  ability  to  analyze  and  predict 
aircraft  behavior  in  those  areas  and  other  areas  that 
are  considered  to  be  at  the  extreme  limits  or  actually 
beyond  the  conventional  flight  envelope. 
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Even  though  such  high  aoa  manuevering  capabilities  exist 
in  some  of  our  current  aircraft,  flight  manuals  often 
restrict  pilots  from  operating  in  this  flight  regime  due 
to  our  inability  to  adequately  control  or  even  understand 
aircraft  behavior  at  such  high  angles  of  attack.  Such 
restrictions,  while  presently  necessary  for  safety 
reasons,  decrease  an  aircraft's  capabilities  in  an  area 
where  an  increase  in  capability  is  vitally  needed  to 
insure  air  superiority  in  air  combat  situations  of  the 
future. 

P.  Overview 

In  chapter  two  the  research  project  will  be 
described.  This  chapter  will  include  a  detailed 
description  of  the  STOL  aircraft,  a  discussion  of  the 
flight  regime  and  types  of  manuevers  to  be  studied  along 
with  a  section  covering  problems  associated  with  high  aoa 
flight  dynamics  analysis.  Chapter  three  will  cover  the 
preliminary  research  taken  including  a  discussion  of  the 
other  methods  considered  and  the  reasons  for  selecting 
bifurcation  and  catastrophe  theory  techniques.  Next,  the 
theoretical  concepts  will  be  developed  and  then  discussed 
as  they  apply  to  the  development  of  the  procedure.  In 
chapter  six  the  specific  methods  and  numerical  techniques 
used  in  the  computational  application  of  the  theroretical 
procedure  will  be  outlined.  A  summary  of  the  results  of 
the  study  will  be  given  next  followed  by  a  discussion  of 
those  results. 
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Conclusions  from  the  study  accompanied  by  recommendations 

of  how  to  best  utilize  the  findings  of  the  study  will 

\ 

comprise  chapter  nine  and  will  conclude  the  text  portion 
of  the  document.  Information  deemed  overly  detailed  for 
the  smooth  flow  of  the  text  will  be  placed  in  the 
appendix. 
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CHAPTER  II 
PROJECT  DESCRIPTION 

_A^_ Description  of  the  Aircraft  and  Data 
The  F-15  STOL  Demonstrator,  (F-15S),  with  three 
major  exceptions,  is  nearly  identical  in  appearance  and 
aerodynamic  performance  to  the  current  inventory  F-15B 
model.  Those  three  exceptions  are  the  addition  of  a  two- 
dimensional  thrust  vector intf /reversing  nozzle  system,  the 
implementation  of  a  close-coupled-canard  and  the 
integration  of  flight  and  propulsion  controls  by  using  an 
aggressive  computer  aided  stability  system  (CASS).  The 
new  thrust  vectoring  nozzle  system  at  the  time  of  this 
study  has  not  been  integrated  into  the  aircraft,  and  it 
will  not  be  included  in  this  study. 

The  close-coupled-canard  is  located  just  forward  of 
the  wing.  This  creates  very  strong  vortex  coupling 
between  the  canard  and  the  wing  and  explains  the  title 
given  to  this  configuration.  The  vortex  coupling  in  much 
the  same  way  as  leading  edge  flaps  or  wing  leading  edge 
extensions  enables  the  aircraft  to  achieve  much  higher 
coefficients  of  lift,  typically  30%  to  60%  higher 
depending  on  the  Mach  number  and  a  stalling  angle  of 
attack  that  is  approximately  three  times  higher  than  the 
unmodified  F-15B.  When  the  aircraft  does  stall,  the 
canard  tends  to  stall  first  causing  a  breakup  of  the 
vortex  structure  across  the  wing. 
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This  in  turn  causes  the  wing  to  stall.  At  this  point  a 
pitch-down  moment  results  which  is  in  most  oases  a 
desirable  effect  [8].  The  canard  is  located  on  the 
engine  inlet  shelves  at  an  average  distance  of 
approximately  2  feet  forward  of  the  root  of  the  wing  and 
is  set  at  a  20  degree  dihedral  angle.  It  has  a  swept 
planform  with  a  sweep  back  angle  of  47  degrees  and  a 
symmetrical  cross-section.  The  canard  airfoil  itself  was 
originally  the  horizontal  tail  of  an  F/A-18.  Apart  from 
the  canard,  all  other  shapes  and  dimensions  are  the  same 
as  the  unmodified  F-15.  Further  details  on  the  physical 
characteristics  of  the  aircraft  are  contained  in  Appendex 
A. 

The  F-15S  flight  control  system  is  hydraulically 
powered  and  fully  irreversible.  Each  of  the  control 
surfaces  have  the  capability  to  not  only  move 
symmetrically  but  differentially  as  well.  This 
significantly  enhances  the  aircraft's  flexibility 
allowing  it  to  have  direct  force  manuevering 
possibilities  as  well  as  reconfigured  controls 
capabilities.  Differential  control  deflections  also 
enable  the  aircraft  to  inorease  its  capabilities  in  the 
combat  performance  areas  discussed  in  chapter  one. 
Perhaps  the  most  important  contribution  of  the  fully 
differential  control  surfaoes,  however,  is  to  give  the 
CASS  more  flexibility  in  selecting  the  correct  control 
positions  to  both  maintain  aircraft  stability  and  follow 
the  commands  of  the  pilot. 


The  CASS  is  a  full  state  feedback  type  control  system 
that  takes  as  inputs  all  of  the  state  variables  in 
addition  to  aircraft  altitude,  mach  number,  pilot  inputs 
and  the  current  positions  of  each  of  the  control 
surfaces.  From  all  of  this  information  the  CASS  selects 
the  control  positions  that  will  result  in  the  performance 
desired  by  the  pilot,  and  it  updates  the  control 
positions  40  times  a  second.  Through  the  use  of  full 
state  feedback  concepts  the  CASS  is  designed  to  maintain 
the  system  eigenvalues  within  certain  boundaries 
selected  to  create  the  desired  aircraft  behavior  for  the 
various  situations  expected  to  be  encountered.  This 
means  that  similar  to  stable  aircraft  without  stability 
augmentation  the  eigenvalues  will  have  negative  real 
parts  or  at  least  if  the  real  parts  are  positive  their 
magnitudes  will  be  small.  The  system  works  quite  well. 
In  fact,  without  the  CASS  the  F-15S  would 
catastrophically  diverge  from  straight  and  level,  trimmed 
flight  in  less  than  ten  seconds  when  left  on  its  own. 
Unfortunately,  the  CASS  is  not  able  to  completely 
maintain  stability  under  all  situations.  He  are 
therefore  left  with  the  typical  situation  of  an  aircraft 
that  is  stable  when  subjected  to  a  large  variety  of 
disturbances  from  equilibrium  yet  unable  to  maintain  its 
stability  under  all  possible  conditions.  The  F-16 
aircraft  is  a  good  example  of  this  situation. 
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It  too  is  controlled  in  a  similar  fashion  and  would  be 
unstable  without  its  CASS.  The  F-16  is  highly 
manueverable  and  under  most  conditions  very  stable  to 
fly.  However,  at  very  high  angles  of  attack  during 
extremely  transient  air  combat  manuevering  the  F-16  has 
on  several  occasions  departed  controlled  flight.  This  is 
due  to  the  CASS's  inability  to  deal  with  the  conditions 
created  by  that  type  of  high  aoa  manuevering,  and 
consequently,  restrictions  have  been  placed  on  the  F-16’s 
flight  envelope. 

The  F-15S  is  projected  to  be  capable  of  some  very 
substantial  increases  in  performance  over  its 
conventional  F-15B  counterpart  [9].  The  takeoff  distance 
will  be  27%  less  while  the  landing  distance  is  projected 
to  decrease  by  82%.  These  improvements  are  primarily  due 
to  the  two-dimensional  thrust  vectoring/reversing 
nozzles;  however,  a  portion  of  the  improvement  in  this 
area  is  due  to  the  decreased  approach  and  takeoff  speeds 
resulting  from  the  decreased  wing  loading  afforded  by  the 
canard.  The  roll  rate  is  expected  to  increase  by  24X 
while  the  pitch  rate  is  projected  to  double.  These  two 
major  increases  in  capability  will  pave  the  way  for  large 
increases  in  aircraft  agility  by  allowing  the  aircraft  to 
quickly  change  its  attitude  and  flight  path.  The 
expected  increases  in  coefficient  of  lift  and  stalling 
angle  of  attack  have  already  been  discussed  above. 


The  addition  of  the  canard  is  expected  to  increase  the 
positive  load  factor  capability  to  9g.  The  actual  role 
of  the  canard  in  this  improvement  is  somewhat 
controversial  at  the  present  time,  but  it  can  be  safely 
said  that  the  decreased  wing  loading  provided  by  the 
canard  does  play  a  role  in  this  area  of  increased 
capability.  Also,  by  enabling  the  aircraft  to  be  trimmed 
to  a  more  streamlined  attitude  during  cruise,  the  canard 
is  expected  to  increase  the  cruise  range  by  13%. 
Finally,  the  F-15S  model  should  be  able  to  decrease  its 
deceleration  time  by  up  to  44%  in  certain  situations 
giving  it  more  flexibility  in  an  air  combat  situation. 
The  above  projected  increases  in  combat  manuevering 
capabilities  are  impressive  when  taken  into  consideration 
that  the  current  F-15  air  combat  capabilities  are  perhaps 
the  best  in  the  world. 

The  data  available  on  the  F-15S  for.  this  research 
project  are  wind  tunnel  studies  accomplished  at  Wright- 
Patterson  AFB,  the  aerodynamic  model  developed  by 
McDonnell  Douglas,  the  canard  schedule,  and  the  CASS 
package  also  developed  by  McDonnell  Douglas.  The  wind 
tunnel  studies  include  static  and  forced  oscillation 
data.  Mo  rotary  balance  tests  have  been  done  yet.  The 
wind  tunnel  data  is  in  both  plot  form  and  tabular  form 
and  is  very  extensive.  In  several  cases  the  coefficients 
are  described  as  a  function  of  four  variables. 
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In  order  to  keep  calculation  costs  and  time  down  the  data 
set  was  simplified  by  restricting  the  configurations  to 
be  investigated  and  by  decreasing  the  number  of  variables 
to  be  used  in  the  determination  of  the  coefficients.  The 
aircraft  configuration  was  restricted  to  gear  up,  flaps 
up  and  no  stores.  After  studying  the  data,  it  was 
determined  that  the  data  did  not  vary  appreciably  with 
changes  in  altitude.  It  was  also  determined  that  the 
data  for  Mach  numbers  higher  than  .6  were  not  necessary 
since  a  study  of  high  aoa  flight  dynamics  was  planned. 
Therefore,  a  data  file  for  a  Mach  number  of  .2  and  one 
for  a  Mach  number  of  .6  were  created.  Both  of  the  files 
were  for  an  altitude  of  20,000  feet.  These  two  data 
files  give  the  ability  to  accurately  look  at  both  spin 
type  motions  and  wind-up  turn  type  manuevers.  These 
simplifications  decreased  the  tabular  data  sizes  so  that 
there  were  no  four  dimensional  tables  and  only  a  few 
three  dimensional  tables.  In  order  to  accurately  model 
spin  type  motions  the  need  for  rotary  balance  data  was 
present.  This  need  was  filled  by  obtaining  rotary 


balance 

data 

for 
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F-15  aircraft 

and 

adjusting 
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data 

to 

account  for 
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data 

was 
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that  the  system  of  aircraft  equations  of  motion  did  not 
receive  an  impulsive  input  when  the  rotary  balance 
correction  factors  were  added  in  (aoa=55  deg.)  and 
removed  (aoa=30  deg.). 


A  complete  discussion  of  the  development  and  application 
of  the  rotary  balance  data  is  contained  in  appendix  B. 
The  aerodynamic  model  except  for  the  addition  of  the 
rotary  balance  correction  factors  is  identical  to  the  one 
determined  and  validated  in  simulations  by  the  Air  Force 
and  McDonnell  Douglas.  The  symmetrical  canard  position 
is  independent  of  the  CASS  system  and  inputs  from  the 
pilot;  it  is  a  function  of  aoa  and  Mach  number  only.  The 
differential  canard  setting  is  controlled  by  the  CASS  and 
is  applied  to  the  current  symmetrical  canard  position. 
When  the  CASS  is  activated,  the  the  model  used  in  this 
study  for  the  F-15S  has  six  independent  control 
surfaces.  The  independent  control  surfaces  are:  left  and 
right  rudder,  left  and  right  stabilator,  differential 
canard  setting  and  the  ailerons.  The  CASS  package  is 
implemented  with  no  modifications  to  keep  the  aircraft 
model  as  close  to  the  actual  aircraft  as  possible. 


B.  Discussion  of  the  Flight  Regime  to  be  Studied 
Ever  since  the  beginning  of  air  combat  manuevering, 
air  battles  would  invariably  progress  towards  higher  and 
higher  aoa  manuevering  in  attempts  to  achieve  a  tighter 
turn  radius  than  the  opponent  so  that  a  favorable 
position  for  weapons  delivery  could  be  gained.  Achieving 
higher  turn  rates  requires  in  most  cases  that  an  aircraft 
fly  at  higher  angles  of  attack. 
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Presently,  fighter  aircraft  handling  characteristics 
at  very  high  angles  of  attack  leave  much  to  be  desired. 
The  problems  range  in  severity  from  being  unable  to 
accurately  track  a  target  due  to  wing  rock  and  other  high 
aoa  oscillations  to  total  departure  from  controlled 
flight  resulting  in  loss  of  the  airoraft  and  possibly  the 
pilot.  Some  specific  high  aoa  flight  problems  are 
discussed  below.  At  high  aoa  the  rudder  controls  become 
sluggish  due  to  the  vertical  tail  descending  into  the 
wake  of  the  wing  and  horisontal  tail.  This  reduces 
weathercock  stability  creating  a  situation  ripe  for  a 
nose  slice  type  departure  from  controlled  flight.  The 
low  aspect  ratio,  swept  back  wing  design  of  most  fighters 
make  control  with  the  ailerons  not  only  ineffective  at 
high  aoa  due  to  wing  tip  washout  but  actually  dangerous 
in  some  cases  where  adverse  yaw  characteristics  are 
severe.  Additionally,  swept  wing  aircraft  experience  a 
decrease  in  static  longitudnal  stability  as  a  result  of 
the  wings  stalling  from  the  tips  inward  forcing  the 
center  of  pressure  location  forward  towards  the  center  of 
gravity  location.  At  high  angles  of  attack  the  induced 
drag  due  to  the  slender,  low  aspect  ratio  wings  becomes 
extremely  high  causing  engine  power  to  be  in  many  cases 
totally  inadequate  to  sustain  manuevers  or  even  control 
the  airoraft. 
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At  high  angles  of  attack  the  lateral/directional  and 
longitudnal  modes  experience  a  high  degree  of  coupling 
causing  an  input  in  one  mode  to  excite  very  undesirable 
behavior  in  the  other.  Another  interesting  problem 
peculiar  to  fighter  aircraft  results  from  their  ‘sleek* 
aerodynamic  design  which  tends  to  concentrate  the  mass  in 
the  fuselage.  This  increases  the  aircraft's 
susceptibility  to  roll  coupling  which  is  a  combination  of 
the  inertial  and  kinematic  coupling  that  occur  when 
rolling  an  aircraft.  Roll  coupling  is  aggravated  by 
both  high  angles  of  attack  and  high  roll  rates  [10]  This 
phenomenon  as  well  as  some  of  the  others  discussed  above 
often  can  occur  very  abruptly  in  a  manner  that  can  be 
termed  as  catastrophic.  The  above  problems  in  addition 
to  a  few  others  currently  make  high  aoa  flight  an 
undesirable  flight  regime  to  operate  in. 

When  all  of  the  undesirable  side  effects  of  high  aoa 
flight  discuss  ve  are  considered  as  a  whole,  the  the 
high  aoa  dynamics  of  a  fighter  aircraft  can  be  at  best 
described  as  unpredictable.  This  is  true  regardless  of 
whether  or  not  the  aircraft  has  been  equipped  with  a 
computerised  stability  system.  As  an  example,  consider 
the  F-16  and  the  F-5E  with  the  shark  nose  modification. 
The  F-16  as  mentioned  above  has  demonstrated  a  great  deal 
of  unpredictability  at  high  angles  of  attack  resulting 
in  several  lost  aircraft. 


) 
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But,  the  F-5E/shark  aircraft  has  also  been  shown  to  be 
very  unpredictable  due  to  the  complex  aerodynamics 
brought  into  play  by  the  shark  nose  section.  It  too  has 
departed  controlled  flight  without  warning  resulting  in  a 
loss  of  aircraft.  In  each  of  the  above  cases  a  lack  of 
overall  understanding  of  how  each  of  the  above  high  aoa 
phenomena  might  combine  to  make  up  the  entire  high  aoa 
flight  dynamics  picture  has  been  directly  responsible  for 


the  lost 

aircraft. 

Before  aircraft 

agility 

can 

be 

increased. 

we 

must 

remove  a  great  deal 

of 
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unpredictability 

from 

the  high  aoa 

dynamics 

of 

our 

fighter  aircraft. 

_C,  High  Angle _of_ Attack  Analysis  Problems 
The  primary  problem  with  studying  high  aoa  flight 
dynamics  is  the  nonlinear  characteristics  of  high  aoa 
flight.  The  great  bulk  of  aircraft  dynamics  analysis  has 

been  based  on  several  assumptions  that  become  totally 
unrealistic  when  considering  an  aircraft  flying  at  high 
aoa.  The  most  far  reaching  of  those  assumptions  was  the 
linearisation  of  both  the  equations  of  motion  and  the 
aerodynamics.  The  linearisation  of  the  equations  of 
motion  is  based  on  the  assumption  of  small  perturbations 
about  an  operating  point  through  the  use  of  the 
Poinoare/Liapunov  theorem  (Appendix  C). 
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The  assumption  of  small  amplitude  disturbances  is  a 
questionable  one  in  our  application  but  can  be  made  to 
hold  up  given  a  small  enough  neighborhood  around  the 
point  of  interest.  Unfortunately,  the  eigenstructure  of 
the  system  goes  through  very  large  changes  when  the  aoa 
is  increased  to  the  very  large  values  being  investigated 
forcing  the  eigenvalues  to  cross  the  imaginary  axis  in 
many  cases.  This  partially  invalidates  the 
Poincare/Liapunov  theorem  forcing  the  use  of  a  variation 
known  as  the  Center  Manifold  Theorem  (Chapter  IV)  to  be 
used  in  its  place.  The  Center  Manifold  theorem  requires 
the  use  at  least  in  part  of  the  nonlinear  representation 
of  the  system.  Additionally,  linearisation  of  the 
equations  of  motion  allows  the  use  of  stability 
derivatives  which  necessarily  implies  a  linearisation  of 
the  aerodynamic  data.  This  is  an  extremely  unrealistic 
assumption  when  dealing  with  high  aoa  aerodynamics .  The 
next  major  assumption  made  to  simplify  the  equations  of 
motion  was  to  assume  that  motions  in  the  plane  of 
symmetry  did  not  create  motions  out  of  the  plane  of 
symmetry.  This  decoupled  the  equations  of  motion  into  a 
fourth  order  set  of  longitudnal  equations  and  a  fourth 
order  set  of  lateral/directional  equations  [6].  For  the 
reasons  discussed  above  this  is  also  an  unrealistic 
assumption  when  dealing  with  high  aoa  flight.  Wing  rock 
is  amoving  the  many  examples  contrary  to  that  assumption. 
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The  decoupled,  linearised  set  of  equations  is  not 
completely  without  merit  though.  They  have  and  still 
will  accurately  predict  aircraft  behavior  in  low  aoa 
small  amplitude  disturbance  flight  with  excellent 
accuracy. 

Perhaps  the  greatest  underlying  cause  for  the  high 
aoa  analysis  problems  discussed  above  is  the  highly 
unsteady,  unpredictable,  nonlinear  aerodynamic  flow 
fields  that  surround  an  aircraft  at  high  aoa.  It  is 
those  flow  fields  that  create  the  nonlinearities  in  the 
data  which  is  primarily  responsible  for  the  large 
variations  in  the  location  of  the  system  eigenvalues. 
Consequently,  it  is  the  aerodynamic  flow  fields  that 
ultimately  force  the  governing  equations  for  the  aircraft 
to  be  nonlinear.  Some  specific  aerodynamic  flow  phenomena 
present  when  flying  at  high  aoa  in  the  F-15S  are 
assymetrio  vortex  shedding,  canard -wing  vortex 
interaction  and  large  amounts  of  separated  and  reattached 
flow  fields  [4,8,11].  These  high  aoa  aerodynamic  flow 
phenomenon  oan  have  both  positive  and  negative  effects  on 
an  aircraft's  flight  performance.  Assymetric  vortex 
shedding,  for  instance,  is  currently  a  problem  because  of 
its  unpredictability,  but  should  the  mechanisms  behind 
the  phenomenon  become  better  understood,  the  use  of 
vortex  shedding  could  eventually  be  used  in  a  positive 
way  to  enhance  aircraft  manueverability. 
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The  vortex  interaction  between  the  canard  and  the  wing 
certainly  has  a  positive  effect  on  the  aircraft's 
performance,  but  it  also  complicates  the  problem  by 
contributing  significantly  to  the  nonlinearity  of  the 
situation  as  a  whole.  Although  this  study  will  not  deal 
with  the  fluid  dynamics  associated  with  the  problem,  it 
is  recognised  that  a  better  understanding  of  the  high  aoa 
fluid  dynamics  involved  will  pave  the  way  for  swifter 
progress  in  the  understanding  of  the  flight  dynamics  at 
high  aoa. 

In  addition  to  increasing  the  nonlinearity  and 
complexity  of  the  aircraft's  flight  dynamics,  the  canard 
presents  some  other  problems  to  be  overcome  in  analysing 
and  controlling  the  aircraft's  flight  behavior.  The 
canard  by  greatly  increasing  the  aircrafts  response  to 
control  inputs  allows  the  aircraft  to  generate  extemely 
high  pitch  rates  so  that  the  aircraft  .  can  routinely 
overshoot  its  stalling  angle  of  attack.  The  higher  pitch 
and  roll  rates  also  create  increased  problems  due  to 
inertial  coupling.  As  is  normally  the  case,  the 
increased  capabilities  afforded  by  a  design  modification 
drastically  increases  the  complexity  of  the  problem  for 
the  engineer. 

The  single-most  difficult  problem  in  effectively 
studying  a  system  as  complex  as  this  one  is  the  lack  of  a 
practical  oriteria  from  an  engineering  viewpoint  to 
disoern  global  stability. 
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The  great  bulk  of  the  nonlinear  analysis  techniques  in 
use  today  are  very  local ized,  the  few  that  are  of  a 
global  nature  are  very  difficult  to  apply  to  different  or 
more  complex  systems.  The  specific  advantages  and 
disadvantages  of  the  particular  nonlinear  analysis 
techniques  available  today  will  be  discussed  in  detail  in 
the  next  chapter.  It  can  be  said  at  this  point  however 
that  the  best  way  to  handle  the  analysis  will  include  a 
combination  of  linear  and  nonlinear  methods  under  the 
general  guidance  of  a  more  global  method  that  will 
effectively  combine  several  localized  techniques  into  a 
complete  analysis  package. 
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CHAPTER  III 
PRELIMINARY  RESEARCH 

A.  General  Points 

A  good  deal  of  preliminary  research  was  necessary 
to  not  only  ensure  that  all  options  were  considered  but 
also  to  enable  me  to  become  more  familiar  with  the 
field  of  nonlinear  analysis  in  general  and  specifically 
as  it  applies  to  aircraft  dynamics.  The  types  of 
research  undertaken  were  quite  varied.  A  literature 
search  of  the  nonlinear  analysis  field  including  texts, 
papers,  journals  and  other  periodicals  was 
accomplished;  indepth  discussions  on  the  subject  with 
several  people  involved  in  the  nonlinear  analysis  field 
were  held;  and  actual  tests  of  the  more  promising 
techniques  were  made.  All  of  the  literary  works 
studied  during  the  reseach  that  had  any  bearing  on  the 
selection  of  the  method  of  research  or  the  method 
itself  are  cited  in  the  references.  Those  individuals 
who  gave  of  their  time  and  expertise  to  further  this 
project  are  gratefully  acknowledged  in  the  foreword. 
And  the  specific  techniques  tested  as  well  as  those 
considered  are  briefly  discussed  in  the  remainder  of 
this  chapter. 
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B.  Method  Selection  Considerations 
Most  of  the  considerations  that  were  taken  into 
account  in  the  selection  of  the  method  of  analysis 
have  already  been  alluded  to  in  one  manner  or  another. 
For  continuity  I  will  briefly  summarize  the  selection 
criteria  used.  First  of  all,  the  method  must  be 
applicable  to  the  F-15S;  it  must  be  capable  of 
analyzing  high  order  nonlinear  systems  of  equations; 
and  the  method  must  be  practical.  The  practicality  of 
the  method  will  be  determined  on  the  basis  of  whether 
it  is  useful  from  an  engineering  point  of  view,  and 
whether  the  concept  itself  along  with  the  difficulty  of 
its  application  is  within  the  scope  of  a  master's 
thesis.  A  fairly  important  personal  desire  was  that 
the  concept  be  relatively  unique.  One  of  the  most 
important  considerations  is  that  the  method  must  have 
the  potential  of  yielding  useful  and  meaningful 

results.  Specifically,  the  technique  should  be  general 
enough  to  be  readily  applied  to  other  systems,  and  to 
be  useful  to  the  fighter  aircraft  design  community,  it 
must  have  the  capability  to  handle  complex  systems. 
Finally,  in  order  to  satisfy  one  of  the  primary  goals 
of  this  study,  the  results  must  be  of  a  global  nature. 
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C.  Techniques  Investigated 
A  large  number  of  possibilities  were  studied 
before  deciding  on  bifurcation  and  catastrophe  theory 
techniques  as  the  best  course  of  action  to  take  in 
analyzing  the  high  aoa  flight  dynamics  of  the  F-15S. 
Methods  from  all  of  the  general  areas  of  engineering 
analysis  were  considered.  These  general  areas  included 
methods  from  the  following  fields:  analytical, 

assymptotic,  topological  and  numerical.  Brief 
discussions  of  the  techniques  investigated  are 
presented  below. 

lz..-Uihear.isqtion  Mgtbgds.-L6i.7j.  12, 133 

The  major  problem  with  all  linearization  schemes 
studied  was  the  lack  of  their  ability  to  obtain  global 
information  on  the  system.  The  very  nature  of  the 
linearization  process  indicates  that  the  information  is 
valid  only  in  a  local  area  about  a  specific  operating 
point.  A  major  improvement  over  true  linearization  is 
found  in  quasi-linearization.  In  this  method  the  small 
signal  requirement  is  removed  by  repeatedly  re¬ 
linearizing  the  system  as  the  input  causes  changes  in 
the  nature  of  the  system.  This  allows  an  unlimited 
variation  of  input  signal,  and  by  changing  the  linear 
characteristics  of  the  approximation  the  basic 
nonl inear ites  of  the  system  can  be  modeled. 
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Although,  quasi -linearization  is  more  global  than  a 
pure  linear  approximation,  two  major  drawbacks  cause 
this  method  to  break  down.  One  major  problem  stems 
from  the  fact  that  the  extreme  nonlinearity  of  high  aoa 
flight  dynamics  will  force  quasi -linearization 
techniques  to  become  impractical  at  some  point.  This 
is  due  to  the  inability  to  isolate  operating  areas  that 
are  large  enough  to  be  useful  and  still  small  enough  to 
give  an  accurate  approximation.  Another  severe  problem 
stems  from  the  inability  at  some  operating  points  for 
linear  approximations  to  even  be  valid.  The 
theoretical  basis  for  quasi-linearization  methods  is 
found  in  the  Poinare/Liapunav  theorem  which  is 
discussed  in  appendix  C.  As  a  primary  analysis  tool, 
the  potential  for  useful  results  was  not  present  in  the 
schemes  studied.  However,  this  does  not  mean  that  all 
uses  of  linearization  methods  are  without  merit.  There 
are  a  number  of  instances  in  a  nonlinear  problem  when 
inexpensive,  reliable  results  can  be  achieved  through 
linearization  techniques.  A  thorough  understanding  of 
the  nonlinearities  of  the  problem  is  required,  however, 
before  linearised  approximations  can  be  safely  used. 


This  technique  is  certainly  global  in  nature.  Its 
unusual  approach  is  to  study  the  stability  of  the 
system  in  the  whole  without  looking  at  the  details. 
Without  question,  it  does  enable  one  to  discern  the 
stability  of  a  complex,  nonlinear  system,  but  the 
method  does  have  some  serious  drawbacks  from  a 
practical  standpoint.  One  major  drawback  stems  from 
the  extreme  difficulty  in  nearly  all  cases  except 
relatively  simple,  potential  governed  systems  of 
deriving  a  Liapunov  function.  Another  extreme  drawback 
to  the  method  from  an  engineering  point  of  view  is  the 
tremendous  lack  of  generality  of  the  Liapunov  functions 
themselves.  For  instance,  should  a  Liapunov  function 
be  found  and  an  analysis  developed  for  one  particular 
system,  the  engineer  will  have  made  very  little 
progress  towards  a  solution  of  a  similar  problem  since 
a  different  Liapunov  function  must  be  derived  for  each 
case.  Overall,  the  method  is  very  impractical  for  use 
in  this  particular  study. 

Phase  Plane_IraJeotory  Analysis  [16,17,18,191 

This  technique  has  some  positive  points  to  offer. 


Its  results  are  easily  interpreted  and  they  are 
slightly  more  global  in  their  nature  than  some  of  the 
other  methods  investigated. 
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On  the  other  hand,  due  to  its  specific  topological 
nature,  it  is  difficult  to  apply  phase  plane  type 
analysis  to  large  scale,  multi -variable  systems,  nor 
does  it  lend  itself  readily  to  systems  that  have 
extreme  nonlinearities.  The  method  was  used  with  some 
success  in  the  analysis  of  a  nonlinear  model  of  the 
lateral /directional  mode  of  an  F-lll.  As  a  primary 
analysis  technique,  however,  phase  plane  analysis  falls 
short  of  the  selection  criteria  established,  but  it 
does  have  some  promise  as  a  technique  to  compliment  a 
more  global,  larger  scale  method. 

4. -Assymptotic  Methods  T20. 21. 22. 23. 241 

Assymptotic  methods  have  proven  to  be  very 
effective  in  developing  useful,  accurate  approximations 
to  systems  that  are  analytically  very  complex.  These 
methods  have  been  shown  to  be  effective  in  analysing 
both  nonlinear  and  nonautonomous  systems  of  equations. 
Multiple  Scaling  which  is  a  specific  assymptotic 
approach  initially  showed  much  promise  as  a  technique 
to  actually  solve  the  nonlinear  equations  of  motion 
with  nonlinear  coefficients  in  a  quasi -analytical 
manner  as  opposed  to  numerical  integration.  This 
method  is  unique  in  that  it  has  not  yet  been 
extensively  applied  to  aircraft  dynamics  problems. 


A  generalized  multiple  scaling  technique  was  tested  on 


the  same  F-lll  lateral/directional  model; 
unfortunately,  only  limited  success  was  attained  during 
the  study.  I  still  feel  that  the  method  has  the 
promise  to  make  a  breakthrough  in  this  area,  but  at  the 
present  time  the  determination  of  the  proper  scaling 
factors  is  as  much  of  an  art  as  it  is  a  science  making 
the  method  impractical  given  the  time  constraints  on 
the  completion  of  this  research. 

*  Numerical  Simulations  f 131 

Numerical  simulations  have  and  will  continue  to  be 
used  as  a  major  flight  dynamics  analysis  tool.  Several 
of  the  techniques  discussed  in  this  chapter  are  made 
possible  through  numerical  simulation  techniques,  and 
many  other  methods  are  validated  through  numerical 
simulations.  The  simulations  can  be  quite  accurate 
provided  that  the  aerodynamic  model  correctly  models 
the  actual  aircraft,  and  the  numerical  techniques  are 
sound.  There  are,  however,  some  significant  problems 
with  using  numerical  simulations  as  the  primary 
analysis  tool  in  a  high  aoa  flight  dynamics  study. 
When  high  aoa  flight  is  considered,  the  simulations 
become  valid  only  in  small  areas  around  the  state  space 
point  being  simulated  due  to  very  complex,  nonlinear 
aircraft  behavior. 
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With  this  in  mind  the  amount  of  simulations  required  to 
gain  an  adequate  representation  of  the  aircraft 
dynamics  rises  drastically  making  the  method  quite 
costly.  Another  serious  problem  arises  from  the  fact 
that  as  the  complexity  of  the  flight  dynamics  increases 
it  becomes  very  possible  to  miss  critical  flight 
behavior  when  studying  an  aircraft  through  numerical 
simulations  alone.  On  the  other  hand,  numerical 
simulations  are  very  useful  if  used  in  conjunction  with 
a  more  global  approach.  The  more  global  technique 
should  be  able  to  locate  areas  where  possible  stability 
problems  that  require  further  study  exist.  A  numerical 
simulation  of  the  area  in  question  could  then  be  done 
thereby  making  numerical  simulations  more  cost 
effective  and  more  useful  in  general. 

6.  Describing  Functions  f 13.253 

A  quasi -linear  function  that  describes  the 
transfer  properties  of  the  nonlinearity  is  called  a 
describing  function.  The  major  drawback  to  this 
technique  is  that  the  form  of  the  signal  must  be 
calculated  or  known  in  advance.  The  form  that  is  of 
most  use  in  this  research  area  is  the  sinusoidal  input 
which  is  useful  for  determining  the  existence  and 
locations  of  limit  cycles.  Another  more  subtle 
drawback  of  this  method  is  that  only  the  specific 
questions  asked  of  it  are  answered. 
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This  makes  the  method  far  too  specific  to  be  of  great 
value  in  this  study,  and  it  is  also  indicates  that  the 
results  will  not  be  global  in  nature.  Therefore,  use 
of  the  method  should  be  limited  to  instances  where  it 
is  desired  to  investigate  a  limit  cycle  possibility 
predicted  by  a  more  global  method. 

_IL__SelectiQn  of.  the.  Actual  Technique  to  be  Used 
After  all  of  the  preliminary  research  was 

completed,  a  technique  based  on  bifurcation  and 
catastrophe  theories  was  selected  as  the  primary 
research  technique  [26,27,28,29].  The  following 
discussion  outlines  the  reasons  for  making  that 

particular  choice.  The  method  was  very  practical  from 
several  viewpoints.  First  of  all,  a  great  deal  of 
information  on  the  technique  was  available  to  me.  The 
concept  of  combining  bifurcation  and  catastrophe 
theories  in  the  manner  used  by  this  method  was  devised 
by  Scientific  Systems  Inc.  which  is  a  company  located 
in  Cambridge  Ma.  They  were  willing  to  share  some  of 
their  knowledge  and  experience  gained  during  their 
research.  This  factor  contributed  towards  enabling  me 
to  study  a  realistic  engineering  problem  within  the 
limited  time  frame  allotted  for  my  master's  thesis. 
The  method  is  also  practical  from  an  engineering 
standpoint  since  it  is  readily  applied  to  engineering 
type  situations. 
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This  is  due  to  the  ability  to  present  the  results  in  a 
easy  to  understand  graphical  format  that  relays  a  large 
amount  of  usable  information  in  a  short  amount  of 
time.  The  technique  is  applicable  to  fighter  type 
aircraft,  and  has  already  been  tested  on  an  older 
current  inventory  fighter.  The  concept  is  quite  general 
and  can  even  be  adapted  to  other  uses  in  addition  to 
aircraft  dynamics.  The  concept  is  also  very  unique 
having  been  applied  to  just  a  few  aircraft  and  none 
with  the  degree  of  sophistication  of  the  F-15S 
aircraft.  The  primary  advantage  of  the  method, 

however,  is  its  ability  to  analyze  systems  from  a 
global  viewpoint.  One  is  able  to  get  a  'picture'  of 
the  flight  dynamics  of  an  aircraft  across  a  full 
spectrum  of  control  inputs  in  a  single  gaze.  Another 
important  point  is  that  the  global  nature  of  this 
method  is  not  dimenished  by  the  presence  of 
nonlinearities.  Overall,  the  general  and  global  nature 
of  this  method  make  it  perfect  for  use  as  a  primary 
research  technique  that  will  promote  more  effective  use 
of  the  techniques  discussed  above. 

E.  General  Applications  for  the  Results  of  the  Study 
A  very  important  consequence  of  the  results  of  any 
scientific  study  should  be  a  better  understanding  of 
the  phenomenon  behind  the  driving  desire  for  the 
study. 
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In  this  case  the  phenomenon  is  high  aoa  flight,  and  a 
better  understanding  of  the  combined  effects  of  the 
many  high  aoa  flight  dynamics  phenomena  is  expected  to 
be  gained  through  this  area  of  research.  A  better 
understanding  of  the  nature  and  the  mechanisms  behind 
high  aoa  flight  could  make  it  possible  to  provide  early 
warning  of  impending  occurences  of  the  more  undesirable 
consequences  of  high  aoa  flight.  This  concept  can  be 
likened  to  stall  warning  devices  in  general  aviation 
aircraft.  Another  application  for  the  method  is  model 
structure  determination  and  validation.  This 
particular  aspect  will  not  be  treated  in  this  study. 
The  types  of  application  that  will  be  concentrated  on 
most  in  this  study  will  be  the  development  of 
particular  control  stategies.  One  obvious  control 
strategy  area  is  to  develop  bifurcation  free  control 
laws  for  computerized  flight  control  systems  that  will 
avoid  the  dangerous  areas  by  avoiding  the  nonlinear 
jump  phenomenon  that  normaly  accompanies  bifurcations. 
Another  more  immediately  useful  goal  is  to  determine 
specific  recovery  techniques  for  certain  out  of  control 
situations.  And  finally,  after  the  aircraft  control 
problems  are  better  understood,  the  development  of  air 
combat  tactics  and  techniques  can  be  undertaken  so  that 
full  advantage  is  taken  of  the  newly  realized  high  aoa 
manuevering  capabilities  of  our  fighter  aircraft. 
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CHAPTER  IV 
THEORETICAL  CONCEPTS 

A.  Overview 

Within  this  chapter  the  theoretical  basis  for  the 
method  of  research  will  be  discussed  in  detail.  This 
will  be  done  by  first  introducing  and  discussing  the 
concepts  of  bifurcations  and  catastrophes.  The 
discussion  will  include  a  couple  of  simple  examples  that 
illustrate  the  two  concepts.  Next,  the  four  major 
theorems  that  will  be  used  to  extend  the  concepts  of 
bifurcations  and  catastrophes  to  ‘real  world'  engineering 
applications  will  be  discussed.  For  clarity  and  brevity 
the  abbreviation,  BACTM,  which  denotes  Bifurcation  and 
Catastrophe  Theory  Methodology  will  be  used  to  refer  to 
the  general  techniques  used  in  this  research.  The  name 
is  adopted  from  the  title  of  the  research  conducted  in 
this  same  area  by  Scientific  Systems  Inc. ,  Cambridge  Ma. 

_B,-_  Description  of  the  Underlying  Concents 
As  a  result  of  the  nonlinear  nature  of  high  aoa 
flight,  aircraft  dynamics  exhibit  a  wide  variety  of 
bifurcation  and  catastrophe  related  behavior. 
Bifurcations  are  not  a  new  concept  nor  has  their  study 
been  limited  to  recent  history. 
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They  were  first  studied  by  Poincare  in  the  late  1800' s 
and  later  by  Hopf  as  well  as  others  [38].  Catastrophe 
theory  on  the  other  hand  is  a  relatively  new  topic.  The 
earliest  work  published  on  the  concept  was  a  paper 
written  by  Rene  Thom  in  1968  [30,31].  The  concept  gained 
a  good  deal  of  popularity  initially,  but  later  it  became 
overused  leading  to  attempts  to  apply  catastrophe 
theories  to  a  large  number  of  non -technical  systems  which 
lacked  definite  mathmatical  structure.  Examples  of  such 
subjective  use  of  the  theorem  can  be  found  in  [32].  This 
has  led  to  some  considerable  criticism  of  the  theory 
itself  as  a  result  of  its  inability  to  give  clear 
objective  results  in  those  cases  [33].  This  does  not 
constitute  a  problem  for  the  application  proposed  here 
since  a  definite  physical  and  mathmatical  model  of  the 
system  exists. 

L._  Bifurcations 

As  alluded  to  above,  bifurcations  are  a  nonlinear 
phenomenon.  They  can  occur  in  a  number  of  different 
circumstances  and  in  a  number  of  different  forms. 
Literally  speaking,  a  bifurcation  is  a  point  where 
something  is  divided  into  two  parts  or  branches.  This 
literal  definition  is  very  close  to  the  actual 
engineering  application  of  the  word.  In  our  case,  the 
branches  are  solutions  to  the  equilibrium  surface 
equations  which  will  be  discussed  later. 
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A  simple  bifurcation  occurs  when  a  system  through  the 
variation  of  some  independent  parameter  reaches  a  point 
where  its  behavior  can  assume  one  of  two  different  forms 
for  the  same  set  of  system  parameters.  The  point  where 
the  system  branches  off  into  one  of  the  two  behavior 
modes  is  termed  a  simple  bifurcation  point.  An  important 
consequence  of  such  a  phenomenon  is  that  the  mathmatical 
structure  of  the  system  itself  also  undergoes  a  change  at 
this  point.  This  gives  rise  to  the  concept  of  structural 
stability  which  should  not  be  confused  with  stability  in 
terms  of  the  aircraft  motion  [16]. 

An  example  of  a  simple  bifurcation  is  shown  below, 

f  (y)  y 

a 


a  y 

Fig.  4.1a  Fig. 4. lb 

Figs.  4.1a  and  4.1b  Simple  Bifurcation  Example 
Figure  4.1a  depicts  a  simple  y  and  f(y)  relationship.  The 
shape  of  the  curve  is  due  to  the  nonlinearities  in  the 
system.  Let  f(y)  be  a  potential  governed  function  such 
as  gravity  where  f(y)  is  height  at  some  point,  y. 
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Therefore,  the  system  is  in  equilibrium  whenever 
df/dy=0.  Potential  governed  systems  such  as  gravity  tend 
to  seek  minimum  potential  energy,  and  therefore,  any 
equilibrium  point  that  is  also  a  local  minimum  will  be 
locally  stable.  Point  a  in  figure  4.1a  represents  a 
stable  equilibrium  since  df/dy=0  and  d  f/dy  is  positive 
at  that  point.  Figure  4.1b  shows  the  location  of  points 
of  minimum  potential  energy  in  the  (y,  z)  space  where  z  is 
a  variable  that  runs  perpendicular  to  the  paper  in 
figure  4. la.  When  the  value  of  z  is  increased  beyond 
point  b,  the  structure  of  the  system  changes  form.  There 
are  now  two  local  minimums  in  the  area  under  study  and 
the  total  number  of  possible  equilibrium  positions  in  the 
area  being  studied  has  been  increased  to  three,  two 
stable  and  one  unstable.  The  system,  as  it  attempts  to 
minimize  potential  energy,  must  proceed  along  one  branch 
or  the  other  as  z  is  increased  beyond  point  b.  This  is  a 
simple  bifurcation,  and  point  b  is  a  simple  bifurcation 
point.  Physically,  this  situation  could  represent  a 
marble  rolling  along  a  channel  towards  an  intersection 
created  by  a  smooth  ridge.  The  top  of  the  ridge  is  a 
locus  of  unstable  equilibrium  points,  and  the  bottom  of 
the  two  valleys  created  by  the  ridge  are  loci  of  stable 
erquilibrium  points.  Given  an  exact  set  of  circumstances 
the  marble  could  roll  up  the  divider  and  along  the  ridge; 
however,  the  most  likely  path  would  be  down  one  of  the 
two  valleys  to  each  side  of  the  ridge. 
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This  is  due  to  the  unstable  nature  of  the  equilibrium 
points  along  the  ridge.  The  use  of  terms  from  topology 
is  by  design  since  the  concepts  of  differential  topology 
will  be  of  great  assistance  in  understanding  and 
analyzing  the  results  of  a  study  using  BACTM.  If  we 
assume  that  figure  4. la  is  a  global  picture  of  the 
locations  of  stable  equilibriums  of  the  system  for  given 
values  of  z,  then  it  can  be  said  that  the  system  will 
behave  in  a  predetermined  manner  for  all  y  and  z  up  to  z 
equal  to  b.  Beyond  the  point  where  z=b  the  global  system 
behavior  becomes  dependent  on  toe  values  of  y  and  z  in 
addition  to  several  other  parameters.  This  is  a  common 
consequence  of  nonlinearities  in  a  system.  When  z  is 
to  the  left  of  point  b  in  figure  4.1b,  the  single  branch 
is  globally  stable  while  the  two  stable  branches  located 
to  the  right  of  point  b  are  only  locally  stable.  The 
domains  of  attraction  of  the  two  locally  stable  branches 
are  very  well  defined.  All  points  above  the  unstable 
equilibrium  belong  to  the  domain  of  attraction  of  the 
upper  equilibrium  branch  while  all  points  below  the 
unstable  equilibrium  belong  to  the  domain  of  attraction 
of  the  lower  equilibrium  branch. 

This  simple  example  can  be  successfully  expanded  to 
include  large  scale  higher  order  systems.  An  important 
step  to  make  beforehand,  however,  is  to  properly  define 
the  terms  global  and  local. 
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For  the  purposes  of  this  study,  the  term  global  will 
indicate  that  the  full  spectrum  of  admissible  state  and 
control  variable  values  is  to  be  considered.  The  term 
local  will  be  taken  to  mean  that  the  properties  of  the 
point  or  area  being  considered  are  valid  only  within  some 
region  that  is  smaller  than  the  full  spectrum  of 
admissible  control  and  state  variable  values. 

Normally,  bifurcation  type  behavior  occurs  in  the 
manner  described  above.  On  occasion,  however,  the 
number  of  stable  branches  that  intersect  at  one  point  are 
greater  than  two.  This  situation  will  be  described  as  a 
general  bifurcation.  Another  form  of  bifurcation  that 
will  be  refered  to  in  this  study  is  the  Hopf 
Bifurcation.  In  the  context  of  aircraft  dynamics,  a  Hopf 
Bifurcation  may  occur  when  a  complex  pair  of  eigenvalues 
of  the  linearised  system  cross  the  imaginary  axis 
(frequency  or  Laplace  domain  description).  There  are  two 
possible  outcomes  of  this  situation.  If  the  system  has 
the  right  combination  of  nonlinear  characteristics,  the 


behavior  will 

grow 

into 

a  stable  limit 

cycle,  or 

if  the 

global  system 

is 

well 

described  by 

the  linearised 

representation. 

it 

will 

show  unstable 

growth 

as  is 

predicted  by  linear  analysis.  Further  discussion  of  Hopf 
bifurcations  will  be  given  later  in  this  section. 


In  many  ways  catastrophes  are  similar  to 
bifurcations.  They  both  occur  in  nonlinear  systems  and 
they  both  involve  a  change  in  the  behavior  of  the 
system.  Like  bifurcation  points,  the  points  where 
catastrophes  occur  mark  a  change  in  the  structure  of  the 
overall  system,  and  also  similar  to  bifurcations, 
catastrophes  occur  as  a  result  of  smoothly  varying  some 
input  to  the  system.  Unlike  bifurcations,  however, 
catastrophes  result  in  a  sudden  and  usually  large  change 
in  the  state  of  the  system,  and  in  most  cases 
catastrophes  also  occur  without  warning.  This  explains 
the  name  given  to  the  phenomena.  Catastrophes  are  also 
described  as  nonlinear  jump  phenomena  compared  to  a 
bifurcation  being  a  smooth  branching  of  the  behavior.  An 
example  of  catastrophe  behavior  in  an  aircraft  would  be  a 
situation  where  at  a  certain  point  a  smooth  increase  in 
aileron  deflection  causes  a  sudden  and  large  increase  in 
roll  rate. 

The  following  example  will  serve  to  give  a  better 
understanding  of  catastrophes  and  how  they  relate  to 
bifurcations.  This  example  is  derived  from  the  cusp 
catastrophe  which  is  one  of  Thom's  seven  elementary 
catastrophes  [27,30,31].  Consider  the  nonlinear 
differential  equation  below. 

dx/dt  =  x3  +  Clx  +  C2  (4.1) 
Where  x  denotes  the  state,  and  (C1,C2)  are  control 
variables. 
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Since  the  system  is  a  scalar  one,  the  system  equilibrium 
may  be  represented  as  the  minima  and  maxima  of  the 
potential  function 

§=  .25x4  +  .  5Clx2  +  C2x.  (4.2) 
Equilibrium  points  occur  whenever  dx/dt=0  in  equation 
4.1.  Let  Cl=-3  for  this  discussion.  It  can  be  shown  that 
there  is  only  one  equilibrium  point  for  |C2|>2  while 
three  equilibrium  points  exist  when  (C2|<2.  This  is  a 
consequence  of  the  number  of  real  solutions  to  equation 
4. 1  when  dx/dt=0.  For  instance,  when  C2=0,  the 
equilibrium  solutions  to  equation  4. 1  are  easily  shown  to 
be  (~V3~  ,  O.VT  );  however,  when  |  C2 1  >2  there  is  only  one 
real  valued  equilibrium  solution  to  equation  4.1. 
Differentiating  equation  4. 1  with  respect  to  x  will 
result  in  the  second  derivative  of  the  potential  function 
and  is  shown  below. 

d2  $  /dx2  =  3x2  +  Cl  (4.3) 
Inserting  the  values  of  x  corresponding  to  the 
equilibrium  positions  (-V3’,0,V^3’  )  into  equation  4.3  show 
that  x=0  is  a  local  maximum  and  therefore,  an  unstable 
equilibrium  point,  and  x=±  3  are  both  minima  and 
therefore,  are  stable  equilibrium  points.  This  is  very 
similar  to  the  marble  in  a  bifurcating  channel  example 
discussed  above  in  that  an  unstable  equilibrium  point 
separates  any  two  stable  equilibrium  points. 
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Refering  back  to  equation  4.3  and  setting  d^<£/dx^  to  0, 
it  can  be  seen  that  regardless  of  the  value  of  C2  there 
will  be  inflection  points  in  the  potential  function  at 
x=±  \/-Cl/3  when  Cl=-3  this  corresponds  to  x=  ±  1.  On  the 
other  hand,  the  location  of  the  maxima  and  minima  do  vary 
with  changes  in  C2.  When  C2  equals  +2  or  -  2  and  Cl=-3, 
a  minimum  and  maximum  will  coalesce  at  one  of  the 
inflection  points.  When  this  occurs,  if  the  system  was 
previously  at  the  stable  equilibrium  point  that  is  now 
an  inflection  point,  the  system  will  be  in  unstable 
equilibrium  and  ready  for  a  catastrophe  or  more  simply 
put,  a  jump  in  its  state  to  occur.  A  small  variation  of 
C2  will  cause  the  jump.  This  example  also  shows  that  a 
stable  equilibrium  cannot  coexist  with  an  unstable 
equilibrium.  The  result  of  any  such  coexistance  is  the 
loss  of  stable  equilibrium  at  that  point.  This  situation 
is  described  as  an  unstable  equilibrium  annihilating  a 
stable  one  [27]. 

The  above  discussion  is  demonstrated  graphically  in 
figure  4.2.  Let  XI  and  X2  be  the  locations  of  the  two 
stable  equilibrium  points,  when  they  exist.  The  system 
will  be  a  marble  seeking  a  local  minimum  of  its  potential 
energy.  This  example  varies  from  the  preceding 
bifurcation  example  in  that  the  behavior  of  this  system 
is  studied  in  only  one  variable.  In  Figure  4.2  the 
control  variable,  C2,  is  allowed  to  vary  slowly  from 
C2=-3  to  C2=+3. 


56 


After  each  change  in  C2,  the  system  is  allowed  to  reach 
equilibrium,  dx/dt=0.  As  discussed  above,  only  one 
equilibrium  point  is  possible  when  |C2|>2.  This  is  shown 
in  plots  a  and  g.  In  plot  b  X2  is  created  when  an 
imaginary  pair  of  roots  to  equation  4. 1  converge  on  the 
real  axis.  This  occurs  at  the  inflection  point,  x=-l. 
Plots  c,  d  and  e  show  the  convolution  of  the  potential 
function  as  C2  is  varied  from  -2  to  +2.  Plot  d  depicts 
the  C2=0  example  discussed  above,  and  therefore, 

X2=+V"5"  and  the  unstable  equibrium  point  is  located  at 
x=0.  In  plot  f  the  local  maximum  and  the  local  minimum  at 
XI  coalesce  at  x=+l.  As  C2  is  increased  further,  two  of 
the  three  real  valued  solutions  to  equation  4. 1  break 
away  from  the  real  axis  and  the  equilibrium  point  at  XI 
disappears.  At  that  point  a  jump  in  the  state  of  the 
system  occurs.  Physically,  the  jump  is  demonstrated  by 
the  marble  rolling  down  to  X2.  If  C2  were  now  varied 
from  +3  back  to  -3,  the  behavior  would  qualitatively  be 


the  same. 

but 

the 

jump 

would  occur  as  C2 

is  decreased 

below  -2 

as 

shown 

in 

plot  b  revealing 

a  definite 

hysteresis  effect. 


Figure  4.2  follows  this  page. 


e.  C2=l 


f.  C2=2 


g.  C2=3 


Figure  4,2:  Potential  Function,  $(x,c)  for 
various  values  of  C2  with  Cl=-3 
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Figure  4.3  brings  the  concepts  discussed  above 
together  in  a  three  dimensional  diagram  in  the  (X.C1.C2) 
space.  The  diagram  shows  a  surface  made  up  of  all  of  the 
solution  points  satisfying  equation  4. 1  for  various 
values  of  the  control  parameters.  Cl  and  C2,  when 
dx/dt=0.  The  area  where  three  real  valued  solutions  to 
equation  4. 1  exist  is  depicted  by  the  equilibrium  surface 
folding  over  on  itself.  A  global  representation  of  the 
positions  where  a  jump  in  state  may  occur  is  shown  by  a 
projection  onto  the  control  space  (C1,C2)  of  all  the 
locations  where  inflection  points  coincide  with 
equilibrium  points.  The  projection  of  these  catastrophe 
point  locations  onto  the  (Cl,  C2)  space  gives  a  mapping 
of  the  locations  where  the  number  of  equilibrium 
solutions  change.  In  this  example  the  area  within  the 
cusp  show  where  the  equilibrium  surface  is  triple 
sheeted,  and  the  area  outside  of  the  cusp  show  where  the 
equilibrium  solutions  are  single  valued.  The  boundary 
between  areas  with  different  numbers  of  equilibrum  points 
is  called  a  bifurcation  surface  since  the  equilibrium 
solutions  separate  or  bifurcate  from  each  other  along 
that  surface.  The  combination  of  all  the  bifurcation 
surfaces  in  the  control  space  is  known  as  a  catastrophe 
map  since  it  gives  the  locations  in  the  control  space 
where  catastrophes  are  likely  to  occur.  Catastrophes 
often  occur  when  a  bifurcation  surface  is  crossed  so  that 
the  number  of  local  equilibrium  solutions  decreases. 
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Figure  4.3;  Equilibrium  and  Bifurcation  Surfaces  for  £ (x,c) 
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The  bifurcation  surface  shown  in  figure  4.3  is  found 
by  locating  all  of  the  points  which  are  both  equilibrium 
and  inflection  points.  This  is  accomplished  by  solving 
the  following  pair  of  simultaneous  equations: 


f(x,c)  =  x°  +  Clx  +  C2  =  0 


(4.5) 


df/dx  =  3x*  +  Cl  =  0 


(4.6) 


Eliminating  x  from  the  above  equations  gives  the  equation 
for  the  bifurcation  surface  as 

( 1/27)C13  +  C2((-1/3)C1)3/2  +  1/4C22  =  0  (4.7) 

The  admissible  values  for  the  control  variable  must  be 
restricted  real  values.  This  requires  that  Cl  be 
negative  for  solutions  to  the  bifurcation  equation  to 
exist.  The  variable.  Cl  is  known  as  a  splitting  factor 
since  for  C1>0  no  catastrophes  can  occur  and  for  C1<0 
catastrophes  can  occur.  This  point  is  shown  graphically 
in  figures  4.4a  and  4.4b  which  are  simply  two  dimensional 
cross  sections  of  the  equilibrium  surface  in  figure  4.3 
for  Cl  held  constant  at  a  positive  value  and  a  negative 
value.  Points  a  and  b  in  Figure  4.4b  are  Bifurcation- 


limit  points. 


stable' 


unstable 


stable. 


Fig.  4.4a  C1>0 


Fig.  4.4b  C1<0 


Figures  4.4a  and  4.4b  2-D  Equilibrium  Plots 


The  simple  example  discussed  above  combines  nearly 
all  of  the  concepts  to  be  used  in  this  study.  A  few 
pivotal  theorems  that  will  be  used  to  extend  the  concepts 
revealed  in  this  scaler  example  to  higher  order,  multi¬ 
dimensional  systems  will  be  discussed  next. 

„CJ„£Ls-g.us.s.i.0D .,-fchg  -Main .  Theorems 
Currently,  there  are  four  theorems  that  are  of 
primary  importance  in  applying  the  concepts  of  BACTM  to 
real  world  engineering  systems.  Each  of  those  theorems 
will  be  discussed  briefly  here. 


1.  The  Center  Manifold  Theorem  f27. 34, 35. 361 

The  Center  Manifold  Theorem  is  of  great  use  in 
determining  the  extent  to  which  linear  techniques  can  be 
applied  to  a  nonlinear  system.  It  is,  in  a  sense,  an 
extension  in  the  application  of  the  Poincare/Liapunov 
Stability  theorem  (appendix  C).  Consider  a  nonlinear 
system  described  by  the  nth  order  differential  equation 
dx/dt  =  Kx.fi)  (4.8) 

where  x  is  a  vector  composed  of  n  state  variables  and  fi 
is  a  vector  composed  of  m  control  variables.  Suppose 
that  the  eigenvalues  for  the  system  are  found  by  using 
some  suitable  linearization  technique.  Let  k  be  the 
number  of  eigenvalues  of  the  linearised  system  whose  real 
parts  are  zero,  i  be  the  number  with  positive  real  parts 
and  j  be  the  number  with  negative  real  parts. 
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The  system  can  now  be  separated  into  three  individual 
systems  by  use  of  appropriate  coordinate  transformations 
so  that  the  original  system  can  be  represented  as 

follows: 

du/dt  =  g(u,c)  (4.9) 

dv/dt  =  A( -) v  +  B(-)c  (4.10) 

dw/dt  =  A(  +  )w  +  B(+)c  (4.11) 

where  (-)  and  (+)  indicate  systems  created  from  the 
portions  of  the  original  system  whose  eigenvalues  have 
negative  and  positive  real  parts  respectively.  The 

systems  represented  by  equations  4.10  and  4.11  are  of 
order  j  and  i  respectively  and  are  properly  described  by 
the  linear  approximations.  This  is  a  consequence  of  the 
Poincare/Liapunov  theorem  (Appendix  C).  The  system 
described  by  equation  4.9  is  a  kth  order  system  created 
from  the  eigenvalues  with  zero  real  parts.  The  essential 
nonlinearities  of  the  original  system  are  described 
completely  by  equation  4.9,  and  since  bifurcations  occur 
only  in  nonlinear  systems  the  complete  bifurcational 
behavior  can  be  studied  by  analyzing  the  k  dimensional 
u-system  by  itself.  This  greatly  reduces  the  size  of  the 
system  that  must  be  studied.  In  fact  Arnol’d  [36,37]  has 
shown  that  the  maximum  value  of  k  is  determined  entirely 
by  the  dimension  m  of  the  control  vector.  This  stems 
directly  from  the  concept  of  structural  stability 
mentioned  earlier. 
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Arnol’d  states  that  the  number  of  singularities  that 
occur  in  "general  position"  is  a  function  of  the  number 
of  independent  parameters  in  the  system.  For  instance, 
if  only  one  control  variable  were  available,  a  single 
eigenvalue  could  be  placed  on  the  imaginary  axis  in  a 
structurally  stable  manner.  There  are,  under  some 
circumstances,  cases  where  systems  with  only  one  control 
variable  may  have  more  than  one  eigenvalue  on  the 
imaginary  axis  at  the  same  time;  however,  in  these  cases 
a  small  shift  in  the  parameters  within  the  system  will 
cause  the  additional  singularities  to  disappear.  In  order 
to  maintain  two  eigenvalues  on  the  imaginary  axis  in  a 
structurally  stable  way,  three  independent  controls  must 
be  available.  When  the  computer  aided  stability  system 
(CASS)  is  not  operational,  the  F-15S  will  have  only  three 
independent  control  parameters,  and  therefore,  k  will  be 
two  or  less  in  these  situations.  When  the  CASS  is 
operational,  however,  the  number  of  independent  control 
parameters  is  increased  to  six,  and  k  may  now  take  on  a 


maximum 

value 

of 

three.  Systems 

with  reduced 

orders 

greater 

than 

two 

have  not 

been 

studied 

to  any  great 

extent 

and 

are 

not  well 

understood. 

This 

is  not 

necessarily  a  major  problem  since,  statistically 
speaking,  the  occurances  of  the  proper  positioning  of  the 
controls  to  obtain  singularities  in  "general  position" 
become  increasingly  rare  as  k  increases. 
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This  conjecture  is  stated  without  proof  by  Arnol'd  [37], 
but  it  does  make  considerable  physical  sense.  Given  that 
this  is  true,  the  possibility  of  having  to  deal  with 
lesser  understood  reduced  order  systems  of  three  may  be 
quite  small.  Arnol'd  also  states  that  any  information 
obtained  on  the  u-system  can  be  generalized  to  the 
original  n£h  order  system  by  "suspending"  the  linear 
systems  having  the  eigenvalues  with  non  zero  real  parts 
to  the  u-system.  Systems  with  k=l,2  are  reasonably  well 
understood  and  well  documented,  and  the  next  two 
theorems  by  Thom  and  Hopf  define  the  bifurcational 
behavior  for  these  systems. 


This  theorem  by  Rene  Thom  [30]  provides 
classifications  of  the  bifurcational  behavior  of  all 
gradient  type  finite  dimensional  systems  where  m  <  6  (a 
gradient  type  system  is  one  in  which  a  potential  function 
exists).  The  assumption  of  a  potential  governed  system 
means  that  the  Jacobian  matrix  describing  the  linearized 
system  will  be  symmetrical,  and  therefore  have  only  real 
eigenvalues  [27].  For  systems  with  m<6,  there  can  be  no 


more  than  two  zero  eigenvalues.  Therefore,  through  the 
Center  Manifold  Theorem  the  complete  bifurcational 
behavior  of  these  systems  can  be  studied  by  studying 
either  a  one  dimensional  or  two  dimensional  system. 
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The  exact  form  of  the  one  or  two  dimensional  system 
depends  on  the  nature  of  the  higher  order  terms  in  the 
Taylor  Series  expansion  of  g(u,c)  (ie.  how  many  of  the 
higher  order  terms  are  properly  neglected).  Thom  has 
labeled  the  general  form  of  the  catastrophes  with  k=l  as 
cuspoids  and  the  ones  with  k=2  as  umbilics.  The  example 
given  in  the  previous  section  falls  in  the  general  group 
of  cuspoids  since  its  reduced  state  space  dimension,  k, 
is  one.  The  control  space  dimension  of  the  example  was 
two  which  further  classifies  the  example  as  a  member  of 
the  family  of  generic  catastrophes  known  as  cusps.  The 
elementary  catastrophe  models  developed  by  Thom  can  be  of 
great  use  in  understanding  the  behavior  of  nonlinear 
systems.  Through  the  use  of  the  Center  Manifold  Theorem 
and  Thom's  classification  theorems  both  model 
identification  and  prediction  of  aircraft  behavior  is 
possible. 

The  above  discussion  has  been  restricted  to  reduced 
order  systems  made  up  entirely  of  real  eigenvalues.  The 
next  theory  deals  with  the  bifurcational  behavior  of 
reduced  order  systems  with  a  complex  pair  of  eigenvalues. 

3.  Hoof  Bifurcation  Theorem  T27.381 

Systems  that  are  of  the  non  gradient  type, 
occasionally  have  equilibrium  solutions  which  consist  of 
closed  orbits  (limit  cycles). 
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Hopf  has  shown  that  a  family  of  limit  cycles  can 
bifurcate  from  an  equilibrium  point  when  a  complex  pair 
of  eigenvalues  of  the  linearized  system  cross  the 
imaginary  axis  with  non  2ero  speed.  In  the  study  of  Hopf 
bifurcations  the  maximum  number  of  eigenvalues  with  zero 
real  parts  is  extended  to  four  in  the  m<6  aircraft 
situation.  This  is  still  considered  to  be  a  system  with 
k<2  since  complex  pairs  tend  to  move  together  in  the  same 
manner  as  the  single  real  valued  eigenvalues  discussed  in 
the  previous  section  [36,37].  The  cases  where  two  complex 
pair  of  eigenvalues  or  a  complex  pair  and  a  single  real 
valued  eigenvalue  cross  the  imaginary  axis  are  not  well 
understood.  Consequently,  in  most  circumstances,  the 
number  of  independent  control  surfaces  used  in  the  study 
of  Hopf  Bifurcations  are  limited  to  two  to  prevent  the 
possibility  of  having  more  than  a  single  complex  pair  of 
eigenvalues  cross  the  imaginary  axis  at  once. 
Additionally,  Hopf  Bifurcations  can  only  occur  if  the 
remaining  eigenvalues  with  nonzero  real  parts  are  well  to 
the  left  of  the  imaginary  axis.  This  is  not  overly 
restrictive,  since  most  aircraft  dynamics  studies  are 
concerned  with  the  points  where  stability  is  lost  for  the 


first  time. 
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4-. _ Global  Implicit  Function  Theorem  F 39 1 

The  global  implicit  function  theorem  was  developed 
by  Palias  and  has  been  successfully  used  to  determine  the 
uniqueness  of  equilibrium  solutions  in  nonlinear 
networks.  The  applications  for  this  theorem  in  nonlinear 
aircraft  analysis  are  straight  forward.  From  the 
discussions  above,  it  is  obvious  that  it  is  very 
important  to  know  the  cases  where  an  aircraft's 
equilibrium  solutions  are  not  single  valued.  Consider 
the  equation  =  0.  Palais'  theorem  states  that  two 
conditions  must  be  satisfied  for  solutions  of  the  above 
equation  to  be  unique.  First,  the  determinant  of  the 
system  jacobian  matrix,  F=  must  be  nonvanishing 
for  all  (x.fi)  (ie.  no  inflexion  points  coincident  with 
equilibrium  points),  and  second,  the  growth  condition 
described  below  must  be  satisfied. 

I  If (x,  c)  I  l-too  as  |  ix|  |  — *  09 

The  importance  of  the  first  condition  is  clear  from  the 
above  discussion  of  bifurcations,  and  if  satisfied,  it 
implies  the  absence  of  bifucations  and  therefore  the 
absence  of  catastrophes  of  the  types  studied  by  Thom. 
The  second  case  addresses  the  situations  where  Hopf 
bifurcations  may  be  present.  In  general,  a  system  that  is 
properly  described  by  a  linearized  representation  should 
always  satisfy  the  Global  Implicit  Function  Theorem  in 
the  valid  range  of  the  linear  approximation. 
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These  theorems  are  important  for  a  couple  of 
reasons.  First,  the  theoretical  development  of  the 
procedure  used  in  this  study  which  is  discussed  in  the 
next  chapter  is  based  on  the  concepts  of  bifurcations  and 
catastrophes,  and  those  concepts  in  turn  are  applied  in 
the  procedure  through  the  use  of  these  theorems.  More 
importantly,  however,  they  provide  an  excellent  basis 
towards  expanding  and  improving  on  the  methods  that  will 
be  brought  out  in  this  study. 
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CHAPTER  V 

THEORETICAL  DEVELOPMENT 

A.  General  Overview  of  the  Method 
The  basic  method  to  be  used  in  this  research  was 
originally  developed  by  Scientific  Systems  Inc.  for  use 
by  the  Navy  [27,28,29].  The  method,  known  as  BACTM 
(Bifurcation  Analysis  and  Catastrophe  Theory 
Methodology),  was  developed  primarily  for  use  on  an 
aerodynamic  data  set  representing  the  McDonnell  Douglas 
F-4  but  was  also  used  on  several  other  more  simple 
aircraft  models  during  the  early  stages  of  its 
development.  The  remainder  of  this  section  will  cover 
the  general  approach  that  will  be  taken  during  this 
study,  and  then  the  following  sections  will  discuss 
each  of  the  major  steps  in  detail.  References  will  be 
made  to  the  theories  and  examples  discussed  in  chapter 
4. 

As  in  any  type  of  scientific  research,  the  first 
major  step  taken  during  this  study  was  the  formulation 
of  the  problem  in  more  specific  terms.  This  included 
the  identification  of  the  state  and  control  variables, 
the  development  of  the  aerodynamic  and  mathmatical 
models,  and  the  selection  of  reasonable  assumptions  to 
simplify  the  problem.  Several  different  sets  of 
equations  of  motion  were  considered. 
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The  ones  selected  as  best  suited  for  the  study  included 
a  fifth  order  and  two  eighth  order  sets  of  equations. 
The  fifth  order  set  uses  as  its  state  variables:  p 
(roll  rate),  q  (pitch  rate),  r  (yaw  rate),  alpha  (angle 
of  attack)  and  beta  (sideslip  angle).  One  eighth  order 
set  employs  the  wind  axis  variables:  p,  q,  r,  alpha, 
beta,  theta  (pitch  angle),  phi  (bank  angle)  and  V 
(velocity),  while  the  other  one  uses  the  body  axis 
variables:  u  (x-axis  velocity),  v  (y-axis  velocity),  w 
(z-axis  velocity),  p,  q,  r,  theta  and  phi.  See  Appendix 
D  for  the  actual  F-15S  aircraft  models  used.  The  fifth 
order  set  of  equations  was  derived  by  holding  V 
constant  and  assuming  gravity  to  be  zero.  By  letting 
gravity  be  zero  theta,  and  phi  decouple  from  the 
remaining  equations  and  thus,  have  no  influence  on  the 
aircraft  dynamics.  This  set  of  equations  is  most 
accurate  in  the  low  angle  of  attack,  trimmed  flight 
range  and  was  used  in  the  preliminary  testing  of  the 
method.  Mehra  and  Carroll  [28]  found  that  setting 
gravity  to  zero  did  not  appreciably  effect  flight 
trajectories  in  the  low  aoa  flight  range  but  did  create 
quantitative  errors  in  tiie  high  aoa  flight  range.  This 
problem  was  especially  prevalent  during  simulations  of 
spin  motions.  Surprisingly,  however,  the  qualitative 
results  of  simulations  of  spin  motions  with  and  without 
the  gravity  terms  did  not  differ  greatly  [27]. 
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The  eighth  order,  wind  axis,  set  of  equations  was  used 
for  increased  accuracy  in  the  higher  aoa  range  during 
BACTM  computations,  and  the  other  eighth  order,  body 
axis  set,  of  equations  was  used  in  the  numerical 
simulations.  A  number  of  additional  variables  were  also 
calculated  during  the  simulations,  and  these  are  listed 
and  discussed  in  Appendix  E.  For  further  discussion  on 
the  mathmatical  models  used  see  Appendix  D.  The 
significant  simplifications  made  in  the  development  of 
the  aircraft  models  are  listed  below: 

a.  Altitude  was  held  constant  during  the  equili¬ 
brium  surface  calculataions. 

b.  The  effects  of  spinning  rotors  were  neglected. 

c.  The  thrust  vector  was  assumed  to  be  in  line 
with  the  x  body  axis. 

d.  Thrust  was  assumed  constant  during  individual 
computer  runs  and  was  set  to  zero  during  spin  type 
motions . 

e.  The  configuration  was  restricted  to  gear  up, 
flaps  up  and  no  stores. 

f.  The  CG  location  was  assumed  to  be  constant 
during  individual  computer  runs. 

g.  The  data  set  was  reduced  to  include  only  the 
data  for  20,000  feet  and  Mach  0.2  or  Mach  0.6. 

The  above  simpl  if icat.ions  were  not  overly  restrictive 
in  that  they  do  not  restrict  the  critical 
nonlinearities  of  the  problem. 
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The  above  restrictions  are  also  quite  compatible  with 
conditions  normally  encountered  in  high  aoa  fight. 

After  the  problem  was  properly  specified  and  the 
aircraft  mode.ls  selected,  equilibrium  surfaces  were 
generated  for  various  control  inputs.  This  involves 
setting  the  time  derivatives  to  zero  in  the  equations 
of  motion  and  finding  the  set  of  solutions  across  a 
spectrum  of  admissible  control  inputs.  At  each  solution 
point  of  the  equilibrium  surface  a  local  stability 
analysis  is  accomplished.  The  movement  and  location  of 
the  eigenvalues  when  correlated  with  the  equilibrium 
surface  can  reveal  a  wealth  of  information  on  the 
nonlinear  nature  of  the  aircraft  model  in  the  area 
under  study.  Some  of  the  information  that  can  be 
gained  is  the  location  of  possible  limit  cycles, 
hysteresis  effects,  domains  of  attraction  of  various 
stable  branches,  locations  of  bifurcations  and  Jump 
phenomenon,  areas  that  are  assymptotically  stable  or 
unstable,  and  steady  state  flight  trajectories.  From 
the  location  of  the  individual  bifurcation  points  which 
were  the  inflection  points  in  the  simple  potential 
example  in  Chapter  IV,  bifurcation  surfaces  are  found. 
These  surfaces  are  the  projections  of  the  equilibrium 
surface  bifurcation  points  onto  a  two  dimensional 
control  space. 


These  bifurcation  points  are  also  the  points  where  the 
eigenvalues  of  the  linearised  system  cross  the 
imaginary  axis,  and  it  is  at  these  points  that  the 
system  either  jumps  to  a  new  equilibrium  position, 
grows  to  a  stable  limit  cycle  or  simply  becomes 
unstable  [15].  The  complete  set  of  all  bufurcation 
surfaces  in  the  control  space  is  referred  to  as  a 
catastrophe  map,  and  it  depicts  the  locations  of 
possible  jump  behavior.  Once  the  results  have  been 
obtained,  predictions  of  the  global  aircraft  behavior 
can  be  made;  stability  criteria  can  be  formulated;  and 
some  possible  bifurcation  free  control  laws  can  be 
developed.  Each  of  these  determinations  can  be  tested 
by  numerical  simulation.  Areas  requiring  further  study 
as  determined  by  the  numerical  simulation  can  then  be 
studied  by  repeating  the  above  process. 

B.  Equilibrium  Surface  Calculations 
As  mentioned  above,  a  large  amount  of  information 
can  be  gained  from  knowing  the  complete  set  of 
equilibrium  surfaces  within  the  admissible  range  of 
state  and  control  variables.  An  equilibrium  surface 
will  not  reveal  information  on  aircraft  behavior  during 
nonlinear  jumps,  but  it  will  give  the  general  direction 
of  the  jump,  the  before  and  after  steady  state 

trajectories,  and  the  location  of  the  jumps. 


Hysteresis  effects  can  be  found  by  simply  reversing  the 
movement  along  the  equilibrium  surface  (in  a  manner 
similar  to  the  minimum  potential  example  shown  in 
figure  4.2).  Just  as  in  the  examples  in  Chapter  IV, 
the  domains  of  attraction  of  the  stable  branches  are 
determined  by  the  locations  of  the  unstable  equilibrium 
branches  (see  figure  4.4b).  Locating  the  domains  of 
attraction  for  a  stable  limit  cycle  is  considerably 
more  complicated,  but  information  towards  that  end  can 
be  gained  by  superimposing  maximum  state  variable 
oscillation  amplitudes  which  are  obtained  by  numerical 
simulation  onto  the  corresponding  points  of  an 
equilibrium  surface.  It  also  should  be  noted  that 
there  can  be  several  different  equilibrium  surfaces 
within  the  admissible  range  of  the  state  and  control 
variables.  This  accounts  for  the  existance  of  separate 
spin  and  trimmed  flight  equilibrium  surfaces.  Strictly 
speaking,  surfaces  considered  separate  in  this 
application  may  not  actually  be  separate  since  they  can 
meet  at  points  that  are  outside  of  the  admissible 
control  and  state  range. 

The  calculation  of  the  equilibrium  surface 
requires  the  solution  of  a  nth  order  nonlinear 
algebriac  equation  for  various  control  inputs.  This  can 
be  done  either  by  solving  for  the  control  positions 
given  a  set  of  state  variables  or  vice  versa. 


The  former  of  the  two  is,  in  most  cases,  easier  due  to 
the  smaller  set  of  equations  and  the  fact  that  many 
aircraft  models  employ  linear  control  variable 
coefficients.  In  the  case  of  the  F-15S,  however,  the 
control  variable  vector  is  not  only  large  (6X1)  but 
also  has  nonlinear  coefficients.  It  is  also  more 
systematic  from  a  research  standpoint  to  select  a 
control  input  as  the  independent  parameter  and 
determine  the  resulting  state  variables  as  the  control 
input  is  varied,  and  therefore,  that  will  be  the 
approach  taken  here.  The  variation  of  a  single 
parameter  to  generate  a  continuous  locus  of  solutions 
from  a  single  known  solution  point  is  known  as  a 
continuation  technique.  This  type  of  technique  is  well 
suited  for  the  computation  of  both  the  equilibrium  and 
the  bifurcation  surfaces. 

Using  state  space  vector  notation,  the  equations 
of  motion  for  the  aircraft  can  be  represented  as 

dx/dt  =  f (x,£) .  (5.1) 

An  equilibrium  surface  by  definition  requires  that 
dx/dt  =  0,  and  therefore,  the  equation  of  an 

equilibrium  surface  is 

1(25. fi)  =  0  (5.2) 

where  x  is  the  state  vector  and  q  is  the  control 
vector.  The  F-15S  aircraft  model  used  in  this  study 
will  require  only  three  control  inputs  although  the 
control  vector  may  be  as  large  as  (6  X  1). 
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If  the  computer  aided  stability  system  (CASS)  is 
operational  then  the  three  control  inputs  will  be 
longitudnal  stick  deflection,  lateral  stick  deflection 
and  rudder  pedal  deflection  as  measured  from  within  the 
cockpit.  The  CASS  will  take  those  ’pit1  inputs,  the 
current  control  positions  and  the  current  aircraft 
state  and  select  the  actual  control  surface 
deflections.  The  number  of  independent  control  inputs 
available  to  the  CASS  in  this  study  is  six.  If  the  CASS 
is  not  operational,  the  number  of  independent  control 
inputs  will  be  decreased  to  three.  In  order  to 
accomplish  this,  the  canard  will  be  scheduled 
symmetrically  as  a  function  of  aoa  and  Mach  number,  the 
stabilator  and  rudders  will  be  operated  symmetrically 
while  the  aileron  surface  movements  remain  unchanged. 
Therefore,  when  computations  are  accomplished  without 
the  CASS,  the  control  inputs  used  in  the  equilibrium 
surface  calculations  and  the  actual  system  control 
inputs  will  be  equivalent.  During  the  equilibrium 
surface  calulations,  two  controls  will  be  held  constant 
while  the  third  is  varied.  This  makes  the  control 
vector,  c,  in  equation  5.2  a  scalar,  and  therefore, 
continuation  techniques  as  described  above  are  easily 
applied.  The  final  form  of  the  equilibrium  surface 
equation  is  then 


f  ( 25.  c)=0 


(5.3) 
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The  first  step  in  an  equilibrium  surface 
calculation  is  to  find  a  single  equilibrium  point 
satisfying  equation  5.3  from  which  to  continue  the 
remaining  solutions  from.  This  can  be  done  by  any 
method  which  presents  the  least  difficulty.  If  a  'trim' 
equilibrium  surface  is  desired,  then  setting  the 
oontrol  inputs  as  well  as  some  of  the  state  variables 
to  zero  will  make  the  analytical  solution  of  an 
equilibrium  point  feasible.  For  a  more  complicated 

equilibrium  surfaoe  .  such  .  as  one  Involving  spin  type 

/ 

motions,  a  steady  state  spin  condition  that  'is  already 
known  to  exist  or  that  is  found  through  simulation  can 
be  used  as  a  ,  smarting  .  point.  Finally,  once  an 
equilibrium  surface  is  generated,  any  point  on  that 
surface  can  be  Used  as  a  starting  point  for  another 
surface.  This  makes  the  method  self  contained  from  the 

•  ,  s 

selection  of  the  initial  equilibrium  point. 

There  are  two  primary  olasses  of  ’continuation 
methods  [28].  The  first  class  is  known  as  continuation- 
by-differentiation  while  the  second  class  is  a  group  of 
iterative  continuation  techniques.  The  continuation -by- 
differentiation  technique  solves  the  equation 

djf/dc  =  F(dx/dc>  +  r  0  (5.4) 
where  as  before  x  is  the  state  vector,  c  is  a  soalar 
control  parameter  and  £  is  £(x#c),  the  right-hand-side 
of  equation  5.3. 


F  is  the  system  Jacobian  matrix 

F  =  (  dfU.c)/**)  (5.5) 

which  can  be  thought  of  as  a  linearized  representation 
of  the  local  system.  Equation  5.4  indicates  that  c  is 
to  be  varied  in  such  a  way  so  that  equation  5.3  will 
always  be  true.  The  major  problem  with  this  method  is 
the  necessity  of  calculating  both  the  Jacobian,  F,  and 
the  vector  d£/bc  at  every  point.  This  becomes  quite 
time  consuming  when  n  is  large.  On  the  other  hand, 
iterative  approachs  use  a  locally  convergent,  iterative 
method  such  as  a  Newton-Raphson  technique  [40,41]  which 
requires  only  the  computation  of  the  Jacobian  matrix  at 
each  point.  The  solution  is  converged  upon  by 
iteratively  solving  an  equation  of  the  following  form: 

2£(new)  =  x(old)  -  F  ^x(old)  (5.6) 

The  convergence  of  the  Newton-Raphson  method  is 
quadratic  providing  that  the  original  estimate  is  near 
the  root  [40].  When  the  estimate  is  not  'near'  the 
root,  the  method  will  usually  fail.  Iterative  schemes 
such  as  the  Newton-Raphson  technique  fall  into  a 
general  class  of  numerical  techniques  known  as 
predictor-corrector  methods  where  the  Newton-Raphson 
formula  shown  above  is  used  as  the  ‘corrector’  step. 
The  use  of  the  Newton-Raphson  technique  alone  is 
therefore  unsuited  for  computation  of  the  equilibrium 


surfaces. 
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Current  continuation  methods  use  a  combination  of  the 
two  classes  discussed  above. 

As  can  be  expected,  the  continuation  schemes 
discussed  above  tend  to  break  down  in  the  neighborhood 
of  bifurcation  points.  It  is  at  these  points  where  the 
slope  of  the  equilibrium  surface  with  respect  to  the 
control  variable  becomes  infinite  or  two  or  more 
equilibrium  surfaces  intersect.  In  a  multi -dimensional 
system  these  occurences  are  indicated  by  the  system 
Jacobian  matrix  becoming  singular.  If  the  Jacobian  is 
non-singular,  then  the  point  is  a  regular  point  with 
£(x,  c),  and  all  of  its  first  partial  derivatives  being 
continuous.  The  implicit  function  theorem  (Chapter  IV) 
ensures  the  existence  of  a  unique  regular  solution  to 
equation  5.3  through  any  regular  point.  Which  in  turn 
ensures  the  success  of  either  class  of  continuation 
techniques  at  those  points.  Conversely  though,  both 
types  of  algorithms  break  down  in  the  neighborhood  of 
bifurcation  points  since  F  becomes  non-invertable  at 
those  points.  Before  proceeding  further  it  will  be  of 
use  to  define  the  three  different  types  of  bifurcation 
points  found  during  the  continuation  process.  The 
first  type  is  called  a  limit  point;  the  second  type 
will  be  called  a  simple  bifurcation  point;  and  the 
third  type  is  known  as  a  general  bifurcation  point.  An 
example  of  each  of  these  is  shown  in  figure  5.1  on  the 
next  page. 


c 

Figure  5.1  Bifurcation  and  Limit  Points 


In  this  study,  a  predictor-corrector  continuation 
method  developed  by  Kubicek  in  1976  [28]  will  be 
employed.  The  Kubicek  algorithm  uses  a  combination  of 
the  continuation-by-differentiation  and  iterative 
methods  discussed  above.  In  dealing  with  the  singular 
Jacobian  matrix,  study  of  the  work  in  this  area  by 
Keller  [42]  gives  much  insight  to  the  nature  of  the 
problem.  Keller  states  that  the  types  of  bifurcation 
points  that  occur  in  the  equilibrium  surface  are 
directly  related  to  the  co-rank  (n-rank)  of  the  system 
Jacobian.  At  regular  points  their  are  no  bifurcations 
and  thus,  the  co-rank=0.  At  limit  points  as  well  as 
bifurcation  points  the  system  Jacobian  is  singular,  but 
the  co-rank  of  the  Jacobian  varies  between  each  of 
these  cases.  Limit  points  have  a  co-rank=l,  and  simple 
bifurcations  have  a  co-rank=2,  etc. 
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During  equilibrium  surface  calculations,  the  state 
vector  is  augmented  with  the  control  vector  selected 
to  be  the  continuation  parameter  so  that  the  state 
vector  becomes  a  (n+1  X  1)  vector  where  n  is  the 
dimension  of  the  system.  The  Jacobian  of  the  augmented 
system  then  becomes  an  (n  X  n+1)  matrix.  At  a  limit 
point  the  Jacobian  of  both  the  original  and  the 
augmented  system  become  singular,  but  within  the 
augmented  system  Jacobian  there  must  be  at  least  one 
nonsingular  (n  X  n)  matrix.  This  nonsingular  matrix  is 
found  by  eliminating  the  most  singular  column  of  the 
augmented  system  Jacobian  resulting  in  a  matrix  of  the 
form  below 

pV4xl . dfl/dxk-l,ifl/9xkn . dV3xn+ll 


Gk  = 


j/dx1 


*fn/dxl'  '  *  •  *^fn/dxk-l’afn/axk+l’  *  *  *  *  *fn/dxn+l 


where 


n+1 


=  c. 


There  are  n+1  possible  (n  X  n)  G, 


matrices  to  choose  from.  The  nonsingular  matrix  G,  is 

XV 

then  substituted  in  place  of  the  system  Jacobian  in  the 
continuation  scheme. 
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In  essence,  the  roles  of  the  control  parameter  and  the 
kth  state  variable  are  now  reversed  so  that  the  state 
variable  corresponding  to  the  column  removed  to  create 
now  becomes  the  continuation  parameter,  and  the 
control  parameter,  c,  now  acts  as  a  state  variable. 
This  situation  remains  intact  until  the  limit  point  is 
traversed  and  the  kth  column  is  no  longer  the  most 
singular  column.  Where  bifurcation  points  are 
encountered,  as  is  often  the  case  in  bifurcation 
surface  calculations  (due  to  the  compression  of  a  three 
dimensional  representation  onto  a  two  dimensional 
surface),  the  state  variable  vector  must  be  augmented 
with  yet  another  control  parameter.  This  process  can 
be  continued  until  all  possibilities  are  exhausted. 
From  the  discussion  in  Chapter  IV,  Arnold  [36,37]  has 
shown  that  there  will  always  be  more  control  variables 
available  than  bifurcation  points  to  contend  with  when 
dealing  with  the  number  of  independent  control 
variables  present  in  this  case  (up  to  six).  Kubicek's 
algorithm  uses  a  lower  order  continuation-by¬ 

differentiation  model  (a  variable  order  Adam-Bashforth 
method)  as  a  predictor  step  and  the  Newton -Raphson 
method  as  the  iterative  corrector  step.  Ferziger  [40] 
does  an  excellent  job  in  explaining  these  two 

techniques,  and  therefore,  they  will  not  be  discussed 
here  in  any  detail. 
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Theoretically,  the  concept  is  complete  at  this 
point.  Unfortunately,  there  are  some  practical  and 
numerical  problems  that  still  must  be  dealt  with.  First 
of  all,  the  actual  application  of  such  a  procedure 
suggested  above  would  be  very  impractical  when  higher 
order  bifurcations  occur  due  to  the  increased  size  of 
the  Jacobian  that  must  be  calculated  at  each  step,  and 
second,  some  numerical  difficulties  still  exist  in  the 
neighborhood  of  bifurcation  and  limit  points  due  to  the 
behavior  of  x(c)  near  those  points.  This  state 
variable  behavior  problem  comes  into  play  in  the 
calculation  of  ?)!/  Sc  in  the  basic  continuation  scheme 
derived  from  equation  5.4  and  shown  below 

dx/dc  =  -F“1(  di/  d  c)  (5.7) 
Kubicek's  method  handles  this  problem  by  introducing  an 
arc  length  parameter,  s,  to  evaluate  the  dependence  of 
the  state  variable,  x»  on  the  control  parameter,  c, 
[28].  By  letting  (X»c)  become  (x(s),c(s))  and 
combining  this  with  the  concept  developed  by  Keller 
which  was  described  above,  the  continuation  scheme 
shown  in  equation  5.7  becomes 

dXj/ds  =  -Gk-1  (dl/ d  xk) (dxk/ds)  j=k,  lljin+1  (5.8) 
where  the  kth  column  is  the  one  replaced  by  the  n+1 
column  containing  the  control  parameter  derivatives.  A 
standard  Euclidean  arc  length  relationship  is  selected 
for  s.  This  is  shown  in  equation  5.9  on  the  next  page. 
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(dx1/ds)2  +  _  +(dxn/ds)2  +  (dc/ds)2  =  0  (5.9) 

By  using  equation  5.8  in  equation  5.9,  dx^/ds  can  be 

found.  The  remaining  n  parametric  derivatives  are  found 

by  back,  substitution  into  equation  5.8.  The  sign 

2 

ambiguity  presented  by  the  solution  of  (dx^/ds)  in 
equation  5.9  is  resolved  by  the  orientation  of  the  arc 
length  parameter,  s,  along  the  solution  branch.  This 
is  determined  at  the  beginning  of  the  algorithm  when 
the  initial  continuation  parameter  direction  is 
selected. 

The  problem  of  dealing  with  higher  order 
bifurcation  points  has  not  specifically  been  addressed 
at  this  point,  but  the  continuation  method  discussed  so 
far  is  quite  adequate  for  the  computation  of 
equilibrium  surfaces  with  limit  points.  There  are 
several  possible  means  to  handle  the  computation  of  the 
equilibrium  surfaces  at  simple  and  general  bifurcation 
points.  One  such  method  has  already  been  alluded  to. 
The  state  vector  could  be  augumented  by  additional 
control  parameters  so  that  a  non-singular  Jacobian 
matrix  can  be  found  in  the  same  manner  described 
above.  This  method  would  work,  but  it  would  also  be 
quite  costly  since  a  larger  dimension  system  would  have 


to  be  dealt  with. 


A  more  comprehensive  arc  length  normalization  scheme 

such  as  that  suggested  by  Keller  [28,42]  would  be 

better  suited  for  the  computation  of  equilibrium  branch 

bifurcation  points.  Unfortunately,  that  type  of  arc 

length  normalization  would  be  overly  complicated  when 

dealing  with  limit  points  alone,  and  experience  has 

shown  that  in  the  computation  of  equilibrium  surfaces 

limit  points  occur  far  more  frequently  than  bifurcation 

points.  The  approach  that  will  be  taken  during  this 

study  will  be  to  adapt  the  Kubicek  method  with  minor 

improvements  (those  improvements  will  be  discussed 

further  in  the  next  chapter).  Bifurcation  points, 

where  two  or  more  equilibrium  branches  intersect,  can 

be  stepped  over  by  simply  using  a  predictor-corrector 

technique  that  will  converge  to  the  other  side  of  the 

bifurcation  point.  In  the  case  of  a  simple  bifurcation 

point  the  procedure  for  locating  the  second  branch  is 

discussed  below.  Study  of  the  two  branches  eminating 

from  the  simple  bifurcation  point  have  shown  that  they 

both  are  tangent  at  the  simple  bifurcation  point  to  the 

plane  described  by  the  two  eigenvectors  associated  with 

T 

the  two  zero  eigenvalues  of  G  G,  which  is  a  square  (n+1 
X  n+1)  matrix  of  co-rank=2  (G  is  the  augmented 
Jacobian)  [28].  The  plane  can  be  determined  from  those 
two  eigenvectors  and  a  search  for  all  points  satisfying 
the  equilibrium  equation,  5.3  can  be  accomplished. 


86 


Once  again,  it  should  be  remembered  that  when  the 
computer  aided  stability  system  is  not  operational,  the 
number  of  zero  eigenvalues  is  limited  to  two,  and 
therefore,  the  equilibrium  surface  bifurcational 
behavior  is  restricted  to  simple  bifurcations  and  limit 
points.  When  the  computer  aided  system  is  operational, 
general  bifurcations  are  possible  but  are  limited  to 
the  cases  where  three  branches  intersect.  As  was 
discussed  in  Chapter  IV,  the  situations  where  the 
reduced  order  system  size,  k  (k=co-rank) ,  is  greater 
than  two  are  not  well  understood.  Locating  the 
additional  equilibrium  branches  eminating  from  a 
general  bifurcation  point,  however,  can  be  accomplished 
in  the  same  manner  described  above  except  that  it  will 
be  progressively  more  difficult  due  to  the  increased 
dimensions  of  the  space  that  must  be  searched  for  the 
points  satisfying  equation  5.3.  For  example,  one 
branch  can  always  be  found  by  the  predictor-corrector 
technique  discussed  above.  The  second  branch  in  the 
case  of  a  simple  bifurcation  point  is  found  by  making  a 
one-dimensional  search  in  circular  coordinates  about 
the  bifurcation  point  in  the  plane  described  by  the  two 
critical  eigenvectors.  In  the  case  of  three  branches 
eminating  from  a  bifurcation  point,  the  eigenvectors 
corresponding  to  the  three  critical  eigenvalues  form  a 
three  dimensional  space. 


6? 


Therefore,  a  two-dimensional  search  in  spherical 
coordinates  must  be  made  about  the  bifurcation  point  to 
find  the  solutions  to  equation  5.3. 

A  considerable  amount  of  detail  has  been  given  to 


the  calculation 

of 

equilibrium 

surfaces 

and 

the 

development  of 

the 

continuation 

algorithm 

used 

in 

BACTM.  This  is  due  to  the  pivotal  importance  of  the 
equilibrium  surfaces  to  the  method  in  general,  and  it 
is  also  due  to  the  universal  applicability  of  the 
continuation  method  to  other  areas  in  this  study  as 
well  as  possible  extensions  to  this  study. 
Specifically,  the  continuation  method  of  Kubicek  will 
be  modified  slightly  and  applied  to  the  bifurcation 
surface  calculations  in  section  D. 

C.  Local .  Stability-Calculations 

'it 

At  each  solution  point  of  the  equilibrium  surface 
a  local  stability  analysis  is  accomplished.  Since  the 
system  *is  an  autonomous  system  (the  coefficients  are 
not  time  dependent),  linearization  of  the  system  at 
each  point  is  valid  and  will  reflect  the  local 
stability  behavior  of  each  state  variable  provided  that 
the  restrictions  of  the  Poincare/Liapunov  theorem 
(Appendix  C)  are  followed.  More  importantly,  though, 
the  local  stability  analysis  will  show  the  movement  and 
type  of  eigenvalues  present  in  the  linearised  system. 
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The  importance  of  this  information  has  already  been 
discussed  and  will  only  be  mentioned  briefly  here.  The 
type  of  eigenvalues  (stable  or  unstable)  determine 
whether  an  equilibrium  branch  is  locally  stable  or 
unstable.  Eigenvalues  crossing  the  imaginary  axis 
constitute  bifurcations  and  the  number  of  those 
critical  eigenvalues  determine  the  size  of  the  reduced 
order,  nonlinear  system  as  defined  by  the  center 
manifold  theorem.  The  reduced  order  nonlinear  system 
can  then  be  studied  through  the  elementary  catastrophe 
theorem  or  the  Hopf  bifurcation  theorem.  The  location 
of  equilibrium  points  that  have  zero  eigenvalues  also 
determine  the  bifurcation  surfaces.  The  local  linear 
analysis  can  also  provide  a  direct  link  to  phase  plane 
analysis  as  well  as  other  topological  studies  of  the 
system.  And  finally,  the  local  stability  analysis  gives 
insight  to  the  possible  applications  and  limitations  of 
linearization  and  quasi-linearization  schemes. 

Mathmatically,  the  local  stability  analysis  is 
accomplished  by  finding  the  characteristic  roots  of  the 
Jacobian  matrix,  F,  for  each  equilibrium  solution 
point.  The  eigenvalues  are  then  characterized  by  their 
location  in  the  imaginary  plane.  Different 
designations  are  given  to  each  of  the  different  types 
of  local  instabilities. 
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In  this  manner  the  type  and  number  of  eigenvalues  that 
cross  the  imaginary  axis  as  the  system  goes  unstable 
locally  can  be  determined,  and  from  this  information  an 
understanding  of  the  systems'  bifurcational  behavior 
can  be  gained. 


D.  Bifurcation  Surface  Calculations 
As  has  been  thoroughly  discussed  earlier,  the 
bifurcation  surfaces  are  intimately  related  to  the 
equilibrium  surfaces  in  that  they  are  merely  a 
projection  of  the  limit  and  bifurcation  points  of  the 
equilibrium  surfaces  onto  a  two  dimensional  control 
space.  It  has  also  been  brought  out  that  the  location 
of  bifurcation  points  can  be  found  in  a  number  of 
ways.  The  method  that  is  most  consistent  with  the 
equilibrium  surface  calculation  method  discussed  above 
is  to  locate  the  points  that  are  both  equilibrium 
points  and  have  a  singular  Jacobian  matrix.  Thus,  the 
equation  of  a  bifurcation  surface  is  [28] 

g(y»ci>=  [f (y, ci ) ,det(F) ]T  =  0  (5.10) 

•  where  y  =  [x,c.]T  and  (c. ,  c.)  are  two  separate  control 

J  l  j 

inputs.  Any  two  control  inputs  can  be  chosen  to  form 
the  control  space  and  the  third  is  held  constant.  One 
of  the  selected  controls  is  allowed  to  vary  as  a  state 
variable  while  the  other  is  varied  as  the  continuation 
parameter. 
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The  solution  of  equation  5.10  yields  a  curve  c.(c.),  c, 

1  J  K 

fixed,  in  the  control  space  that  is  a  locus  of  the 

bifurcation  point  projections.  The  complete  set  of  such 
bifurcation  curves  is  known  as  a  catastrophe  map.  The 
system  described  by  equation  5.10  is  simply  the 
equilibrium  system  augmented  by  the  constraint, 
det(F)=0,  and  therefore  has  the  dimension  of  n+1  where 
n  is  the  order  of  the  system.  This  accounts  for  the 
addition  of  a  control  parameter  as  an  extra  dependent 
variable.  From  this  point  the  solution  of  this  system 
of  equations  proceeds  identically  to  that  of  the 

equilibrium  surface  equations  with  three  exceptions. 
First  of  all,  the  dimension  of  the  system  is  increased 
by  one,  making  the  calculation  of  the  surface  more 
costly.  Secondly,  the  nature  of  compressing  a  three 
dimensional  locus  of  points  into  a  two  dimensional 
representation  dramatically  increases  the  possibilities 
of  intersections  occuring  between  the  surfaces,  which 
as  discussed  above,  greatly  increases  the  complexity  of 
the  continuation  process.  The  final  exception  is 
caused  by  the  constraint,  det(F)=0.  Due  to  the 
dimension  of  the  system,  it  is  impractical  to  expand 
the  determinant  of  F  analytically;  therefore,  numerical 
differentiation  of  a  system  with  a  singular  or  near 
singular  matrix  is  required  which  is  not  only  tricky 
but  could  be  quite  costly  as  well. 
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E,  Analysis  of  the  Results 
Due  to  their  global  nature,  the  results  of  an 
analysis  using  BACTM  can  first  be  studied  from  a 
qualitative  viewpoint  with  the  purpose  of  better 
understanding  the  nonlinear  phenomena  in  the  area  under 
study.  Then,  more  quantified  results  can  be  obtained 
in  the  areas  that  show  unusual  or  counter-intuitive 
behavior.  This  approach  is  somewhat  unique  in  that 

many  approaches  reverse  the  above  order  of 
investigation.  Some  examples  of  BACTM  results  from  two 
much  simpler  aircraft  models  than  the  F-15S  are 
discussed  below.  Their  coefficients  are  linear  and 
their  data  bases  are  quite  small  when  compared  to  the 
that  of  the  F-15S.  The  results  are  therefore  much 
clearer  and  more  suited  to  a  general  discussion  of  the 
analysis  techniques  used  in  this  study. 

Ll— Bali.  Coupling  Equilibrium  Surface  Example 

One  of  the  better  known  nonlinear  jump  phenomena 
is  roll  coupling  which  was  first  studied  by  William 
Phillips  in  1948.  Analytical  approaches  to  studying 
roll  coupling  are  fairly  complex  and  cumbersome  [43], 
and  as  will  be  shown  below,  the  criteria  obtained 
through  those  traditional  techniques  are  not 
necessarily  accurate  [27].  Roll  coupling  is  a 
phenomenon  that  occurs  primarily  in  modern  fighter  type 
aircraft. 
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This  problem  stems  directly  from  the  low  moment  of 
inertia  about  the  x  axis  that  is  characteristic  of 
fighter  aircraft.  It  is  both  a  kinematic  coupling 
(exchange  of  angle  of  attack  and  sideslip  during  roll) 
and  inertial  coupling  phenomenon  that  occurs  during 
high  roll  rates  and  is  aggravated  by  high  angles  of 
attack.  In  order  to  combat  this  problem  many  fighter 
type  aircraft  have  restrictions  that  limit  the  maximum 
stick  travel  and  g  loading  during  continuous  aileron 
rolls  [10]. 

Figures  5.2,  5.3,  5.4  and  5.5  show  four  different 
equilibrium  surfaces  for  aircrafts  A  and  B,  the  two 
fighter  type  discussed  above.  The  figures  were 
generated  by  Scientific  Systems  Inc.  during  the  early 
development  of  BACTM.  Figure  5.2  shows  an  aileron 
deflection  equilibrium  surface  for  aircraft  A  at  an 
elevator  deflection  angle  of  -0.5  degrees  and  rudder 
deflection  angle  of  zero.  As  should  be  expected  with 
this  simple  model  the  progression  of  roll  rate,  p,  with 
aileron  input  is  very  linear  throughout  the  aileron 
range,  and  the  local  stability  analysis  shows  that  each 
point  is  stable.  This  is  indicated  by  the  S  appearing 
above  the  curve:  Figure  5.3  again  shows  aircraft  A  at 
the  same  conditions  as  Figure  5.2  except  now  the 
elevator  input  has  been  changed  to  0.0  degrees. 
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At  this  point  the  nonlinearities  of  the  aircraft  model 
are  evident  and  are  indicated  by  the  limit  points 
present  in  the  curve.  Figure  5.3  clearly  shows  the 
presence  of  a  jump  in  roll  rate  at  DA=±8  degrees  as 
aileron  input  is  increased  in  either  direction. 

Although  the  jump  at  DA=±8  degrees  is  not  tremendous, 
the  situation  that  has  developed  upon  crossing  that 
*  jump  has  the  potential  to  become  disasterous.  Most 

aircraft  are  protected  against  this  occurence  by  a 
flight  manual  restriction  requiring  that  a  load  factor 
of  one  g  be  maintained  during  multiple  aileron  rolls. 
Unfortunately  however,  the  difference  in  aircraft 
‘feel’  between  one  g  (Fig.  5.2)  and  zero  g  (Fig.  5.4) 
during  continuous  aileron  rolls  is  not  always  obvious, 
yet  the  difference  in  aircraft  behavior  can  be 
catastrophic.  As  an  example,  refer  once  again  to  fig. 
5.3,  suppose  that  aileron  deflection  is  decreased  from 
0  to  -30  degrees  and  back  to  0  again.  During  that 
control  sequence,  the  roll  rate  was  increased  to  about 
225  deg/sec,  and  a  jump  of  about  60  deg/sec  occured  at 
DA=-8.  As  DA  is  decreased  back  to  zero  the  roll  rate 
will  only  decrease  to  about  130  deg/sec.  In  an  attempt 
to  stop  the  roll  rate  the  pilot  increases  aileron 
deflection  further  until  at  about  15  degrees  of  aileron 
deflection,  the  roll  rate  undergoes  a  catastrophic  jump 
to  -225  deg/sec. 
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Closer  inspection  of  the' s ituation  that  has  developed 
shows  that  the  aircraft  cannot  return  to  its  original 
starting  position  of  zero  roll  rate  by  using  aileron 
alone.  This  is  due  to  the  'zero  roll  rate'  equilibrium 
branch  being  embedded  inside  of  the  hysteretic  loop 
shown  in  Figure  5.3.  It  should  be  noted  that  the  other 
states  are  also  undergoing  similar  nonlinear  jumps  at 
the  jump  points  shown  in  these  figures.  A  small 
elevator  change  to  -0.5  degrees  deflection  will  enable 
the  aircraft  state  to  be  returned  to  the  starting 
point.  The  much  used  roll  coupling  criteria,  developed 
by  Phillips,  gives  a  critical  aileron  deflection  value 
of  +15  degrees  for  this  situation  which  corresponds  to 
the  larger  catastrophic  jump  from  +225  to  -225  deg/sec. 
The  BACTM  study,  on  the  other  hand,  reveals  a  more 
useful  critical  aileron  value  of  +8  degrees.  Numerical 
simulations  have  verified  the  BACTM  selected  value  of 
+8  degrees  as  the  critical  aileron  deflection  angle 
[27].  Figure  5.4  show  the  second  aircraft,  B, 
undergoing  an  aileron  change  during  a  pitch-up 
manuever.  The  curve  has  distinct  linear  and  nonlinear 
regions.  The  limits  of  these  regions  are  marked  by  the 
jumps  at  +34  degrees.  Since  the  catastrophic  jump 
occurs  at  the  extreme  limits  of  control  deflection  and 
at  a  load  factor  well  above  one,  the  nonlinear  jump 
shown  in  Figure  5.4  is  safely  avoided  by  adhering  to 
the  flight  manual  restrictions  listed  above. 


Figure  5.5  again  shows  aircraft  B  in  the  same  flight 
conditions  as  Figure  5.4.  In  this  case,  however,  the 
elevator  deflection  is  held  at  +2  degrees  and  similar 
to  Figure  5.3  the  'zero  roll  rate'  equilibrium  branch 
is  again  embedded  inside  of  a  hysteretic  loop. 

This  example  has  brought  to  light  some  of  the 
information  that  can  be  gained  from  the  study  of 
equilibrium  surfaces.  The  steady-state  trajectory  is 
shown  for  a  varying  aileron  input.  The  nonlinear  jumps 
or  catastrophes  are  clearly  shown,  and  the  effects  of 
hysteresis  are  evident.  Additional  information  is 
gained  by  the  local  stability  analysis  which  is 
represented  by  the  letters  indicating  the  local 
stability  of  the  individual  solution  points  along  the 
equilibrium  curve.  As  is  shown  in  the  legend  on  the 
figures  a  U  stands  for  unstable  and  an  S  indicates 
stable.  The  U  designation  also  indicates  that  the 
local  analysis  is  unstable  as  a  result  of  a  single  real 
valued  eigenvalue  in  the  right-half-plane  while  the  L 
designation  indicates  that  the  local  stability  analysis 
is  unstable  due  to  a  complex  pair  of  eigenvalues  in  the 
right-half-plane.  The  boundaries  between  two  types  of 
stability  behaviors  is  indicated  by  a  hash  mark,  !,  on 


the  curve. 
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The  type  of  local  stability  behavior  corresponding  to 
each  section  of  the  curve  is  indicated  by  the  letter 
above  or  below  that  section  of  the  curve.  Refering 
back  to  Figure  5.5,  the  topological  nature  of  the 
stable  and  unstable  equilibrium  branches  is  clearly 
shown.  Between  each  stable  branch  there  is  an  unstable 
branch,  and  the  unstable  branches  clearly  mark  the 
domains  of  attraction  of  the  stable  branches.  The 
local  stability  analysis  in  Figure  5.5  also  shows  that 
the  limit  points  do  indeed  occur  when  a  single 
eigenvalue  crosses  the  imaginary  axis.  Another  point 
brought  to  light  by  the  local  stability  study  is  that 
the  bifurcational  behavior  of  the  system  can  be  studied 
in  terms  of  a  single  variable  which  is  a  direct  result 
of  the  center  manifold  theorem.  Wherever  a  transition 
from  an  S  to  an  L  designation  is  found,  a  Hopf 
bifurcation  is  present  and  a  limit  cycle  may  exist. 
The  lower  case  letters  on  figures  5.4  and  5.5 
correspond  to  points  on  the  roll  rate  catastrophe  map 
discussed  in  the  next  section. 

Figures  5.2  through  5.5  follow  this  page. 
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Aircraft  A  (DR=0,  DE=-.5  deg.) 


S  =  Stable 
U  =  Unstable 
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5.2t  Equilibrium  Roll  Rate  (P)  vs.  Aileron  Angle  (DA) 
Aircraft  A  (DR=0,  DE*0.0) 
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2^_Roll  Coupling  Catastrophe  Map  Example 

Figure  5.6  shows  a  three-dimensional  section  of 
the  roll  rate  equilibrium  surface  discussed  above  for 
DA  (-35,  +35),  DE  (-2,  +2),  and  DR=0.  The  projection 

of  the  limit  points  from  the  equilibrium  surface  onto 
the  elevator  versus  aileron  control  space  will  form  a 
roll  rate  catastrophe  map.  This  is  shown  in  the  lower 
portion  of  Figure  5.6.  The  shape  and  behavior  of  this 
particular  catastrophe  map  is  very  similar  to  one  of 
the  elementary  catastrophes  attributed  to  Thom  known  as 
the  butterfly  catastrophe  for  its  distinct  shape 
[27,30,31].  Figure  5.7  shows  a  more  comprehensive  view 
of  the  catastrophe  map.  The  numbers  represent  the 
number  of  equilibrium  solutions  existing  at  that 
combination  of  controls.  The  horizontal  lines  show  the 
plane  of  figures  5.4  and  5.5.  Nonlinear  jump  behavior 
may  occur  anytime  a  boundary  in  Figure  5.7  is  crossed 
so  that  the  number  of  equilibrium  solutions  decreases. 
High  aoa  catastrophes  can  be  found  by  holding  aileron 
deflection  constant  and  increasing  elevator.  In  this 
case  departure  may  be  expected  at  DA=0  and  DE=9.3. 
This  too  was  verified  by  simulation  as  being  correct 
[27].  Finally,  the  area  where  bifurcations  can  not 
occur  is  shown  as  the  region  where  only  one  equilibrium 
surface  solution  exists.  Control  laws  can  be  developed 
to  remain  in  this  non-bifurcation  region  by  using  a 
combination  of  the  control  surfaces. 


Equilibrium  Surface 


Fig»  5 « 6«  Equilibrium  Surface  and  Catastrophe  Map 
for  Aircraft  B  (DR  »  0) 
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Pig.  5.7»  Catastrophe  Map  in  (DA»  DE.  DR  *  0)  plane  for 
Aircraft  B.  The  numbers  in  each  region  indicate  the 
number  of  equilibrium  solutions.  The  letters  correspond 
to  points  shown  in  Figures  5*^  and  5.5. 
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F.  Numerical  Simulations 

Once  the  results  of  the  study  have  been  analyzed, 
numerical  simulations  are  employed  to  verify  the 
predictions  made  and  the  criteria  developed.  The 
numerical  simulations  use  the  eighth  order  body  axis 
model  as  the  Bet  of  equations  being  actually 
integrated,  but  a  number  of  other  variables  and 
parameters  of  interest  are  calculated  as  well.  Among 
those  additional  values  calculated  are  angle  of  attack, 
sideslip  angle  and  total  velocity  so  that  the 
bifurcation  system  can  be  directly  studied.  Some 
additional  values  of  interest  that  are  monitored  during 
the  simulations  are  total  vehicle  kinetic  energy  and 
rate  of  change  of  kinetic  energy.  These  values  have 
shown  some  promise  as  possible  quantitative  indicators 
of  impending  catastrophes  [28].  The  complete  list  of 
the  additional  values  monitored  during  the  numerical 
simulations  and  their  definitions  is  contained  in 
Appendix  E.  The  numerical  simulation  method  has 
provisions  for  allowing  altitude  to  vary  and  for  the 
varying  of  the  controls  to  simulate  different  control 
strategies.  The  numerical  simulation  technique  also 
has  provisions  for  the  calculation  of  phase-plane 
trajectories  as  well  as  the  standard  time  history 
plots. 


103 


CHAPTER  VI 
NUMERICAL  APPROACH 

A.  Overview 

The  iterative  nature  of  this  method  makes  it  a 
perfect  candidate  for  evaluation  by  numerical  methods. 
In  fact  the  only  way  that  any  useful  amount  of 
information  can  be  gained  through  this  method  is  by  the 
repeated  numerical  evaluation  of  the  values  discussed 
in  the  proceeding  chapter.  The  numerical  steps  taken 
in  this  analysis  as  well  as  the  specific  numerical 
techniques  used  to  accomplish  those  steps  are  discussed 
in  this  chapter.  Several  of  the  techniques  used  are 
well  known  and  will  not  be  covered  in  great  detail,  but 
the  manner  in  which  they  are  combined  to  create  the 
complete  method  will  be  covered  in  sufficient  detail  to 
give  an  understanding  of  the  overall  process. 

IL-JPata.  Preparation 

As  has  been  mentioned  above,  the  data  set  for  the 
F-15S  is  extensive  which  made  it  necessary  to  reduce 
the  number  of  independent  variables  contained  in  the 
data.  This  was  accomplished  by  restricting  the 

altitude  to  20, 000  feet  and  the  Mach  number  to  either 
Mach  0.2  or  Mach  0.6.  After  these  restrictions  were 
applied,  the  data  set  was  reduced  to  tables  that  were  a 
function  of  one,  two  and  three  variables. 
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In  order  to  ensure  a  smooth  variation  of  the 
aerodynamic  coefficients  as  the  aircraft  state  changes, 
a  curve  fit  that  is  smooth  at  all  points  is  required. 
Cubic  spline  approximations  were  selected  as  the  best 
approach.  Cubic  splines  represent  the  data  by  placing 
a  third  order  polynomial  between  each  data  point  and 
then  matching  the  first  and  second  derivatives  of  the 
adjoining  curves  at  each  data  point.  This  ensures  the 
smoothness  of  the  partial  derivatives  of  the 
aerodynamic  coefficients  with  respect  to  the  state  and 
control  variables.  For  the  tables  that  were  a  funtion 
of  two  variables,  bi-cubic  spline  interpolations  were 
used.  This  is  similar  to  the  cubic  spline  method 
except  that  a  cubic  polynomial  in  two  variables  is  fit 
between  the  four  data  points  that  represent  a  section 
of  the  two-dimensional  data  table.  Bi-cubic  splines 
have  all  of  the  properties  of  cubic  splines.  Since 
the  nature  of  the  first  and  second  derivatives  at  the 
end  points  of  the  tables  were  not  precisely  known, 
natural  cubic  splines  were  selected.  A  natural  cubic 
spline  forces  the  second  derivative  to  be  zero  at  the 
end  points  thus  minimizing  the  overall  curvature  of  the 
curve.  The  spline  coefficients  which  represent  the 
coefficients  of  the  cubic  polynomial  between  each  data 
point  were  calculated  separately. 
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Since  the  values  of  the  spline  coefficients  will  remain 
unchanged  so  long  as  the  data  remains  unchanged,  the 
spline  coefficients  only  need  to  be  evaluated  once. 
The  spline  curves  are  then  evaluated  each  time  a 
variable  that  is  a  function  of  the  table  approximated 
by  the  spline  changes.  During  a  dynamic  simulation, 
all  of  the  values  are  evaluated  on  each  call  to  the 
aerodynamic  coefficient  subroutine,  but  during  the 
determination  of  the  partial  derivatives  for  the 
Jacobian  matrix,  only  values  that  are  varied  during 
that  call  to  the  subroutine  need  be  calculated.  This 
is  discussed  in  more  detail  below.  The  cubic  and 
bicubic  spline  curves  are  generated  and  evaluated  by 
International  Mathmatics  and  Science  Library  (IMSL) 
routines.  The  remaining  tables  that  are  functions  of 
three  variables  were  interpolated  by  linear 
interpolation.  The  difficulty  and  expense  of  applying 
a  spline  curve  fit  in  three  variables  out  weighed  the 
possible  curve  continuity  problems  that  could  arise. 
Additionally,  the  number  of  tables  requiring  linear 
interpolation  and  their  overall  effect  as  compared  to 
the  other  tables  on  the  aerodynamic  coefficients  was 
small.  All  of  the  interpolation  methods  maintain  the 
values  at  the  boundary  of  the  table  when  an  independent 
variable  exceeds  the  bounds  of  the  table  (ie.  no 
extrapolation  was  accomplished). 


Ferziger  [40]  discusses  the  origin  and  application  of 
splines  as  well  as  a  possible  computer  algorithm  for 
their  implementation. 


In  this  section,  I  will  briefly  describe  the  basic 
flow  of  the  algorithm  to  give  some  order  to  the 
description  of  the  specific  techniques  that  follow. 
The  first  major  step  after  the  initialization  of  a 
number  of  quantities  necessary  for  the  control  of  the 
routine  and  the  calculation  of  the  equations  of  motion 
is  to  read  in  the  spline  fit  data  and  parameters  for 
the  specific  Mach  number  of  interest  (Mach=0.2  or 
Mach=0.6).  The  system  Jacobian  is  then  calculated  and 
the  Newton -Raphs on  technique  is  used  to  converge  on  the 
initial  point  from  the  given  initial  conditions.  Then 
the  Adams-Bashforth  method  is  used  to  predict  the  next 
point.  The  Newton-Raphson  technique  is  then  used  again 
to  converge  to  the  next  solution  point.  During  the 
Newton-Raphson  process,  the  augmented  system  Jacobian 
is  calculated  and  Gaussian  elimination  with  controlled 
pivoting  is  used  to  select  the  continuation  parameter 
and  solve  for  F~*.  This  step  is  repeated  until 
convergence  occurs.  Once  convergence  is  achieved 
within  a  given  criteria,  the  Adams-Bashforth  method  is 
again  used  to  predict  the  next  point. 


This  process  is  continued  until  a  state  or  control 
parameter  exceeds  the  limits  set  prior  to  the  start  of 
the  program.  Then,  if  desired,  the  routine  returns  to 
the  original  starting  point  and  generates  a  branch  in 
the  opposite  direction.  The  computation  of  a  complete 
branch  is  rarely  required  in  the  case  of  the  F-15S 
model  due  to  the  symmetric  nature  of  its 
lateral/directional  data. 

2  -Jacob  ian.. .  Ca  1  c.ul  afci  on 

The  calculation  of  the  augmented  system  Jacobian 
is  accomplished  by  finding  the  partial  derivative  of  £ 
with  respect  to  each  state  variable  and  the  control 
parameter,  c.  Due  to  the  complexities  of  the  data, 
numerical  techniques  were  used  to  evaluate  the  partial 
derivatives.  This  task  is  accomplished  by  first 
evaluating  the  aerodynamic  coefficients  at  the  initial 
state  and  control  value,  (x.c),  and  then  using  them  in 
the  equations  of  motion  to  find  £(;&,  c).  In  the  same 
manner,  the  value  of  £  is  found  at  four  additional 
points  by  varying  a  state  variable  twice  in  each 
direction  from  its  original  value  and  then  evaluating  £ 
at  each  of  those  points.  A  cubic  spline  curve  fit  is 
then  made  to  each  of  the  resulting  five  point  curves 
(there  is  one  curve  for  each  of  the  n  components  of 
£).  This  passes  a  set  of  n  smooth  curves  through  the 


original  component  values  of  £. 
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Each  curve  shows  the  variation  of  £  as  a  function  of 
the  state  variable  that  was  varied  to  generate  that 
curve.  The  derivatives  *  of  each  of  the  n  curves  with 
respect  to  the  state  variable  varied  to  generate  them 
are  then  easily  found.  This  process  is  repeated  for 
each  of  the  n  state  variables  and  the  one  control 
parameter  resulting  in  the  (n  X  n+1 )  augmented  Jacobian 
matrix  shown  in  chapter  5.  Some  checks  are  made  to 
ensure  that  the  steps  made  to  each  side  of  the  state 
variables  create  a  meaningful  change  in  £  so  that  an 
accurate  partial  derivative  can  be  found.  The  change 
in  £  must  be  large  enough  to  be  representative  but  not 
so  large  as  to  lose  accuracy. 

3.  Corrector  Phase 

As  has  been  mentioned  above,  the  corrector  step 
used  in  this  algorithm  employs  the  Newton-Raphson 
method  shown  below. 

h  =  -Gk_1f  (6.1) 

where  h  isan(nXl)  vector  of  corrector  steps  and 
is  the  Jacobian  matrix  of  the  augmented  system  with  the 
kth  column  removed.  The  kth  column  represents  the  most 
singular  column  in  the  augmented  Jacobian  matrix.  This 
is  selected  by  Gaussian  elimination  with  controlled 
pivoting. 
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This  process  involves  selecting  a  pivot  element  that 
has  the  largest  magnitude  of  all  candidate  elements. 
Then,  all  elements  that  share  the  same  row  or  column 
as  the  pivot  element  are  reduced  to  zero  by  matrix  row 
operations.  The  above  process  is  repeated  until  n  rows 
and  columns  have  only  one  element.  Mo  pivot  element 
may  be  used  more  than  once  in  this  process.  The 
remaining  column  is  the  kth  column  that  is  removed  as 
the  most  singular  one.  If  n  non-zero  elements  cannot  be 
found  in  this  manner,  then  the  rank  of  the  augmented 
Jacobian  is  less  than  n  and  6^  is  non-invertable.  The 
pivoting  is  controlled  in  the  Kubicek  algorithm  by 
selectively  scaling  the  variables  [28].  This  is 
accomplished  by  weighting  each  of  the  n+1  columns  in 
the  augmented  Jacobian  matrix.  The  column  associated 
with  the  variable  selected  as  the  continuation 
parameter  is  given  the  smallest  weighting  factor.  In 
this  way  the  variable  selected  as  the  continuation 
parameter  is  ensured  to  remain  as  the  continuation 
parameter  unless  an  actual  singularity  is  encountered. 
When  G  is  reduced  to  cannon i cal  form  by  the  above 
process,  it  is  easily  inverted,  and  then  equation  6.1 
is  solved  for  the  Newton-Raphson  corrector  step.  The 
Jacobian  for  the  new  values  of  x  found  through  the 
corrector  step  is  then  calculated.  The  above  process  is 
repeated  until  a  solution  point  satisfying  a  preset 
precision  criteria  is  found. 
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In  order  to  increase  efficiency,  an  option  has  been 
included  in  the  . routine  to  avoid  the  recalculation  of 
the  Jacobian  on  each  step  as  the  method  approaches  the 
solution  point.  The  logic  behind  this  option  stems 
from  the  fact  that  the  steps  made  in  the  Newton-Raphson 
method  become  very  small  when  the  solution  point  is 
approached,  and  therefore,  the  Jacobian  will  not  vary 
much  between  those  steps.  Since  the  Jacobian  matrix  is 
totally  a  function  of  the  state  variables  and  the 
controls,  the  Jacobian  is  calculated  only  when  x  varies 
a  predetermined  amount  from  the  point  where  the 
Jacobian  was  last  calculated.  The  Jacobian  is  also 
always  calculated  at  the  beginning  of  the  search  for  a 
new  solution  point.  Not  calculating  the  Jacobian  at 
each  step  tends  to  increase  the  total  number  of 
iterations  required  to  reach  the  solution  point  but  at 
the  same  time  decreases  the  overall  time  to  reach  the 
solution  point  [40].  The  Gaussian  elimination 
technique  described  above  is  well  known.  Further 
information  on  the  subject  can  be  found  in  [40,41]. 

4.  Local  Stability  Analysis 

Once  convergence  to  a  solution  point  has  been 
accomplished,  the  system  Jacobian  and  state  variables 
at  that  point  are  retained.  Through  the  use  of  an  IMSL 
subroutine  the  eigenvalues  and,  if  desired, 
eigenvectors  are  calculated. 
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The  eigenvalues  are  classified  by  the  type  and  number 
of  unstable  eigenvalues.  The  eigenvectors  can  be  of 
use  when  the  rank  of  G  is  less  than  n.  This  occurs  in 
the  higher  order  bifurcation  situations  discussed  in 
chapter  5.  Consequently,  the  eigenvector  calculation 
option  is  usually  put  into  effect  during  the 
bifurcation  surface  calculations. 

5.  Predictor  Step 

The  inverse  of  the  Jacobian  matrix,  G,  is  also 
saved  from  the  last  iteration  of  the  Newton-Raphson 
corrector  method.  It  is  then  used  in  equations  5.8  and 
5.9  to  solve  for  the  n+1  parametric  derivatives, 
dxi/ds.  After  this  information  is  obtained,  the 
variable  corresponding  to  the  kth  column  is  stepped  a 
predetermined  amount.  The  remaining  n  parametric 
derivatives  are  then  used  in  an  Adams -Bashforth  method 
to  predict  the  next  point  on  the  curve.  Adams -Bash forth 
methods  are  implicit,  multi-step  methods  whose  order 
can  be  varied  to  achieve  the  desired  accuracy.  The 
particular  method  used  in  this  approach  is  a  variable 
order  algorithm  with  a  maximum  order  of  four.  The  use 
of  the  Adams -Bashforth  method  over  the  simpler  Adams - 
Moulton,  or  any  other  explicit  method,  ensures  that 
stability  can  be  unconditionally  guaranteed  for 
whatever  order  is  selected  by  the  routine. 


Due  to  the  popularity  of  the  Adams  methods  in  general, 
a  good  deal  of  information  on  them  can  be  found  in  the 
literature  [40,41]. 

The  new  predicted  solution  point,  found  through 
use  of  the  Adams  integration  step,  where  all  n+1  values 
have  been  changed  is  the  next  starting  point  for  the 
Newton -Eaphs on  iterations.  The  complete  process  is 
repeated  until  the  control  surface  deflection  or  one  of 
the  state  variables  exceeds  a  preset  limit.  The  routine 
will  then  generate  another  branch  from  the  starting 
point  in  the  opposite  direction  if  desired.  As  was 
mentioned  above,  this  is  not  usually  necessary  due  to 
the  symmetrical  nature  of  the  F-15S  lateral/directional 
aerodynamic  model. 

The  bifurcation  surface  routine  is  handled  in 
exactly  the  same  manner  except  for  the  addition  of  the 
constraint  that  the  Jacobian  be  singular  at  each 
solution  point.  During  the  computation  of  the 
bifurcation  surfaces,  the  eigenvectors  are  usually 
calculated  as  well  as  the  eigenvalues.  As  was 
discussed  in  Chapter  5,  the  eigenvectors  can  be  used  as 
an  aid  in  the  location  of  the  additional  branches 
emanating  from  simple  and  general  bifurcation  points. 


The  procedure  described  above  will  yield  a  locus 
of  points  that  describe  the  variation  of  the 
equilibrium  solutions  of  each  of  the  state  variables  as 
a  single  control  is  varied.  In  order  to  get  a  more 
complete  picture  of  a  particular  flight  regime  several 
equilibrium  branches  must  be  created  by  varying  the 
other  controls  also.  The  results  of  these  runs  are 
best  understood  when  placed  in  graphical  form.  Three- 
dimensional  plotting  routines  would  present  the  most 
information  in  a  single  plot  but  have  the  draw  backs  of 
being  more  complicated  to  generate  and  more  difficult 
to  obtain  specific  quantitative  information  from. 
Therefore,  the  results  of  this  study  will  be 
represented  in  a  series  of  two-dimensional  plots  that 
are  linked  together  by  indicating  common  points  between 
certain  plots.  Local  stability  information  and  in  some 
cases  maximum  oscillation  amplitudes  are  also  contained 
in  these  plots. 

The  numerical  simulations  are  accomplished  by 
integrating  the  equations  of  motion  through  use  of  a 
fourth  order  Runge-Kutta  method.  A  fourth  order  Runge- 
Kutta  method  is  a  predictor-corrector  method  that 
employs  an  explicit  method  as  the  primary  predictor 
step  followed  by  two  implicit  corrector  steps  and  one 
implicit  predictor  step. 


s  l  ;k% 
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The  method  is  applied  through  the  use  of  an  algorithm 
known  as  DYSYS  which  stands  for  DYnamic  SYstem 
Simulation.  The  numerical  form  used  by  DYSYS  to 
integrate  the  equations  of  motion  is  shown  below: 

Y(tQ  +  dt)  =  YQ  +  1/6  <K0  *  2KX  ♦  2K g  +  K3)  (6.2) 

where  YQ  is  a  vector  of  the  initial  conditions  for  Y, 
dt  is  the  time  step  and  the  K’s  are  defined  below: 

Kq  =  F(tQ,Y0>  *  dt 

K1  =  F(t0  +  dt/2*  Y0  +  1/2<Ko)>  *  dt 

Kg  =  F(t0  +  dt/2,  Yq  +  1/2<K1 > )  *  dt 

K3  =  F(tQ  +  dt,  Yq  +  K2)  *  dt 
where  F  is  a  vector  of  the  first  order  time  derivatives 
of  Y.  The  use  of  implicit  corrector  steps  gives  the 
method  stability,  and  the  order  of  the  method  gives  the 
accuracy  acceptable  for  this  type  of  engineering 
application.  Much  has  been  written  about  Runge-Kutta 
methods,  and  therefore,  I  am  not  going  to  discuss  the 
method  further  here.  More  information  on  the  subject 
can  be  found  in  [40,41]. 

The  complete  package  of  numerical  techniques  used 
in  this  study  has  been  made  as  efficient  and  universal 
as  possible.  Most  of  the  subroutines  are  shared  by  the 
different  programs,  and  little  would  be  required  to 
expand  on  the  methods  currently  developed  or  apply  them 
to  a  different  aircraft  model. 
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CHAPTER  VII 

PRESENTATION  AND  DISCUSSION  OF  THE  RESULTS 
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A.  Outline  of  Work  Accomplished 
A  large  number  of  equilibrium  surface  calculations 
and  numerical  simulations  were  attempted.  As  a  result  of 
the  numerical  and  aircraft  stability  problems 

encountered  during  the  course  of  this  investigation  only 
a  small  portion  of  those  computer  runs  yielded  meaningful 
results.  The  numerical  difficulties  encountered,  among 
other  things,  resulted  in  no  equilibrium  surfaces  being 
successfully  computed  with  thrust  while  the  unstable 
nature  of  the  eighth  order  aircraft  model  reduced  the 
number  of  meaningful  numerical  simulations  to  a  small 
number.  Due  to  the  complexity  of  the  aircraft  model  and 
the  method  itself,  the  equilibrium  surface  calculations 
were  initially,  quite  oostly.  A  number  of  improvements 
on  the  numerical  efficiency  of  the  technique  have  been 
made  during  the  course  of  the  study  so  that  the  computer 
time  expended  during  the  equilibrium  surface  computer 
runs  has  been  reduced  to  a  reasonable  level.  Some 
logistical  problems  as  well  as  certain  time  constraints 
set  on  the  completion  of  this  project  resulted  in  little 
being  accomplished  in  the  area  of  studying  the  airoraft 
with  an  operational  computer  augmented  stability  system. 
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The  work  that  was  accomplished  in  this  area  was  limited 
to  setting  the  system  up  for  operation  on  our  computer 
facilities  and  developing  the  means  for  its  integration 
into  the  existing  algorithms.  Therefore,  all  of  the 

results  reported  here  are  for  the  aircraft  without  the 

% 

computer  aided  stability  system  operational.  Given  the 
high  computer  run  times  experienced  at  the  beginning  of 
this  investigation  as  well  as  the  lack  of  substantial 
bifurcational  behavior  in  the  equilibrium  surfaces,  it 
was  felt  that  calculation  of  bifuracation  surfaces  were 
not  warranted,  and  therefore,  no  bifurcation  surface  runs 
were  attempted.  From  a  computational  point  of  view,  the 
greatest  amount  of  success  was  found  during  the 
calculation  of  the  fifth  order  equilibrium  surfaces. 
Consequently,  seven  fifth  order  equilibrium  surfaces  were 
successfully  computed  while  only  two  eighth  order 
equilibrium  surfaces  were  successfully  computed.  This 
disparity  is  a  direct  result  of  the  numerical 
difficulties  presented  by  the  eighth  order  model.  One  of 
the  eighth  order  equilibrium  surfaces  generated  was  a 
fepin  surface,  and  one  of  the  fifth  order  equilibrium 
surfaoes  showed  distinct  bifurcational  behavior.  These 
two  particular  surfaces  combined  with  one  of  the  other 
fifth  order  equilibrium  surfaces  where  roll  coupling 
effects  were  noted  contain  the  most  interesting  findings 
of  the  study. 
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The  next  section  describes  and  discusses  the  individual 
results  in  detail. 

B.  Equilibrium  Surface  Results 
In  each  of  the  plots  presented  in  this  section  the 
following  symbols  are  used  to  annotate  the  local 

stability  behavior: 

S  =  all  eigenvalues  are  stable 
U  s  a  single  real  valued  eigenvalue  is  unstable 
L  =  a  single  complex  pair  of  eigenvalues  are 
unstable. 

Combinations  of  the  letters  shown  above  simply  indicate 
additional  unstable  roots  in  the  form  indicated.  The 

plotting  routine  used  to  display  the  results  is  a  simple 
* connect  the  dots  with  straight  lines'  type  of  routine, 
and  therefore,  some  roughness  in  the  curves  may  be  noted. 
This  can  easily  be  corrected  by  decreasing  the 

continuation  step  sise  in  the  equilibrium  surface 

algorithm.  In  some  oases  erratio  behavior  of  the  curves 
may  result  from  the  equilibrium  surface  convergence 
criteria.  On  the  other  hand,  however,  some  discontinuous 
behavior  in  the  curves  can  be  attributed  to  the  type  of 
nonlinear  Jump  behavior  sought  after  in  this 
investigation.  The  subtle  differences  in  these  cases 
will  be  pointed  out  below. 
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1.  Fifth  Order  Equilibrium  Surfaces 

During  the  fifth  order  investigation,  all  of  the 
equilibrium  surfaces  were  generated  with  the  velocity 
set  to  500  f/s,  altitude  set  to  10,000  ft  and  thrust  set 
to  sero.  Additionally,  the  data  set  for  Mach=0. 6  was  used 
for  each  of  the  fifth  order  equilibrium  surface 
calculations.  With  the  exception  of  the  first  and 
seventh  fifth  order  cases  the  results  are  presented  in  a 
series  of  six  figures  contained  on  three  pages.  The  first 
page  shows  the  pitch  rate  and  angle  of  attack  plots 
versus  the  appropriate  control  variable;  the  second  page 
shows  the  roll  rate  and  yaw  rate  plots;  and  the  third 
page  shows  the  sideslip  angle  plot  and  the  yaw  rate  plot 
again.  The  repeated  presentation  of  the  yaw  rate  plot 
here  as  well  as  the  repeated  presentation  of  a  few  plots 
during  the  discussion  of  the  eighth  order  equilibrium 
surface  results  is  done  to  enable  easy  comparisons  to  be 
made  between  related  plots.  Plots  that  are  repeated  are 
given  the  same  title  to  avoid  confusion. 

a.  Case  Humber  One  (DAI=DRU=Q.  DEL=+8  to  -30) 

The  first  fifth  order  equilibrium  surface  generated 
was  a  pitch  up  and  pitch  down  manuever  from  DEL=-5 

degrees. 
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During  this  run,  both  the  aileron  and  rudder  deflections 
were  held  to  zero  while  the  stabilator  deflection  angle, 
DEL,  was  varied  from  -5  to  4-8  and  then  again  varied  from  - 
5  to  -30  degrees.  The  results  of  this  computation  are 
shown  in  figures  7. la  and  7. lb.  Due  to  the  symmetrical 
nature  of  the  fifth  order  model  at  relatively  low  angles 
of  attack  the  variables  p,  r  and  beta  were  essentially 
zero  for  the  entire  run  and  therefore  are  not  presented. 
The  plots  in  Figure  7 . 1  show  that  angle  of  attack  and 
pitch  rate  are  very  closely  related,  and  this  point  is 
confirmed  by  inspecting  the  mathmatical  model  used.  The 
model  (Appendix  D)  shows  that  pitch  rate  should  have  the 
dominant  effect  on  the  rate  of  change  of  angle  of  attack 
during  low  angle  of  attack  flight  conditions.  This  is 
confirmed  by  Figures  7.1a  and  7.1b.  The  local  stability 
analysis  shows  a  transition  to  an  *L*  instability  in  the  - 
15  to  -20  degrees  stabilator  deflection  range.  This 
corresponds  to  an  angle  of  attack  just  below  the 
stalling  angle  of  attack  and  quite  possibly  indicates  a 
Hopf  Bifurcation  to  wing  rook  or  some  other  type  of  high 
aoa  oscillatory  phenomenon.  This  could  be  confirmed  by 
numerical  simulation;  however,  a  fifth  order  simulation 
package  was  not  developed.  The  somewhat  discontinuous 
increase  in  both  aoa  and  pitch  rate  at  DEL=-22  could  be  a 
small  oatastrophe  whose  limit  points  were  stepped  over  by 
the  algorithm. 
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In  terms  of  comparison  to  expected  aircraft  behavior,  it 
is  interesting  to  note  that  an  aoa  of  10  degrees  is 
generated  at  0  degrees  stabilator  deflection.  This 
result  was  not  confirmed  one  way  or  the  other  due  to  the 
lack  of  available  data  of  this  type.  Three  of  the  fifth 
order  equilibrium  surface  calculations  reported  below 
share  points  from  this  run  as  their  initial  starting 
point. 

b.  Case  Number  Two  fDRU=0,  DEL=-10.2. DAI =0  to  -30) 

The  second  surface  calculated  was  an  aileron  roll  to 
the  left  at  approximately  24  degrees  angle  of  attack. 
Note  that  the  aileron  deflection  sign  convention  for  the 
F-15S  is  unconventional  with  negative  aileron  deflection 
generating  a  left  rolling  moment.  In  this  set  of  figures 
rudder  deflection  is  held  to  zero  and  stabilator 
deflection  is  held  to  -10.2  degrees.  Figures  7.2a  and 
7.2b  show  that  at  these  conditions  angle  of  attack  and 
pitch  rate  are  not  effected  by  aileron  control  input. 
Figures  7.2c  and  7. 2d  show  that  a  significant  amount  of 
adverse  yaw  is  created  by  the  aileron  input.  This  is 
demonstrated  by  the  right  yaw  rate  that  is  generated  by 
left  aileron  input.  At  this  angle  of  attack  the  loss  of 
aileron  roll  control  authority  is  evident  and  is 
demonstrated  by  a  right  roll  rate  being  generated  by  a 
left  aileron  deflection. 
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This  effect  is  present  to  a  certain  degree  in  all  high 
performance  aircraft  due  to  their  swept back,  low  aspect 
ratio  wings.  The  requirement  for  an  aileron  rudder 
interconnect  that  blends  in  pro-rudder  as  a  function  of 
angle  of  attack  when  aileron  input  is  commanded  is 
verified  in  Figures  7.2c  and  7. 2d.  During  a  similar 
equilibrium  surface  calculation  at  a  low  angle  of  attack, 
a  pro  yaw  situation  was  noted  during  aileron  input,  and 
this  was  confirmed  as  being  one  of  the  aerodynamic  design 
stategies  employed  by  McDonnell  Douglas  to  minimize  the 
effects  of  adverse  yaw.  On  the  next  page  Figures  7.2e 
and  7. 2d  show  that  yaw  rate  and  sideslip  angle  are 
related  in  a  manner  typical  of  uncoordinated  flight 
conditions  in  that  a  positive  yaw  rate  creates  a  negative 
sideslip  angle.  The  local  stability  analysis  showed  all 
points  to  be  stable  during  this  run.  This  is  partially 
substantiated  by  the  linearity  of  the  curves  shown  in 
Figures  7.2a  through  7.2e. 

o.  Case  Number  Three  (DAI=0.  DEL=-19.6.  DRU=0  to  -301 

Figures  7.3a  through  7.3e  show  the  effects  of  a  high 
angle  of  attack  (approximately  30  degrees)  right  rudder 
roll.  In  this  set  of  figures,  aileron  deflection  is  held 
at  zero  and  stabilator  deflection  is  held  to  -19.6. 
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The  initial  point  of  this  run  corresponds  to  a  stable 
point  very  near  the  L/S  stability  boundary  shown  in  case 
number  two.  As  can  be  seen  in  the  figures,  the  solution 
points  throughout  the  rudder  deflection  range  remained 
locally  stable.  Figures  7.3a  and  7.3b  indicate  that  at 
these  conditions  angle  of  attack  and  pitch  rate  are  not 
effected  by  rudder  control  input.  Figures  7.3c  and  7.3d 
reveal  a  very  linear  increase  in  yaw  rate  and  roll  rate 
with  rudder  input.  Figures  7.3e  and  7.3d  on  the  next 
page  show  that  the  rudder  roll  at  these  conditions  is  a 
coordinated  manuever  with  zero  sideslip  angle  throughout 
the  rudder  control  input.  This  is  an  indication  that  the 
rolling  moment  is  generated  by  the  velocity  differential 
seen  by  the  left  versus  the  right  wing  due  to  the  yaw 
rate.  This  is  somewhat  unusual  in  that  sideslip  is 
usually  the  mechanism  responsible  for  the  rolling  moment 
during  a  rudder  roll.  This  set  of  results  confirm  that 
rudder  should  be  the  primary  control  surface  for 
generating  rolling  moment  in  the  angle  of  attack  range 
shown. 

Case  Humber  Four  (DEL=B. 3. DRU=0. DAI=0  to  -30) 

This  equilibrium  surface  is  also  derived  from  the 
first  equilibrium  surface,  and  it  is  a  surfaoe  for  left 
aileron  roll  at  a  negative  angle  of  attack.  In  this  set 
of  figures  rudder  deflection  is  held  to  zero  while  the 
stabilator  defleotion  is  constant  at  5.3  degrees. 
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During  this  aileron  roll  manuever,  the  adverse  yaw 
effects  present  in  case  number  two  are  not  present  but 
evidence  of  roll  coupling  is  present.  This  is  shown  by 
the  increase  in  negative  pitch  rate  and  negative  angle  of 
attack  with  aileron  input  in  Figures  7.4a  and  7.4b.  This 
type  of  behavior  is  caused  by  the  high  roll  rate 
generated  by  the  aileron  deflection  combined  with  the 
aircraft  inertial  axis  being  offset  from  the  axis  that 
the  aircraft  is  attempting  to  roll  about.  Figure  7.4c 
shows  that  the  roll  response  to  the  ailerons  is  good 
throughout  the  aileron  deflection  range.  Figures  7.4d 
and  7.4e  further  confirm  the  existence  of  inertial  roll 
coupling  type  motion  depicting  a  build  up  of  a  positive 
yaw  rate  while  sideslip  angle  remains  slightly  greater 
than  zero  throughout  the  manuever.  This  indicates  that 
the  yaw  rate  component  is  about  the  aircraft  velocity 
vector.  Physically,  the  situation  is  very  similar  to  a 
spin  type  motion.  This  can  best  be  understood  by 
considering  each  of  the  three  rotational  components  as  a 
single  angular  velocity  vector.  At  zero  angle  of  attack 
and  zero  sideslip  angle  (assuming  lxz=0)  the  aircraft 
inertial  axis  is  coincident  with  the  flight  path  of  the 
aircraft.  The  angular  velocity  vector  of  a  roll  about 
the  flight  path  would,  in  this  case,  be  composed  entirely 
of  p,  the  aircraft  x-body  axis  roll  rate. 
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At  angles  of  attack  or  sideslip  angles  other  than  zero, 
however,  an  attempted  roll  about  the  aircraft  x-body  axis 
using  the  ailerons  will  result  in  a  combination  of  rolls 
about  each  of  the  body  axis  with  the  overall  angular 
velocity  vector  still  aligned  with  the  flight  path 
vector.  Rotating  bodies  tend  to  seek  a  mass 

distribution  that  is  symmetric  with  respect  to  its  axis 
of  rotation.  Therefore,  the  heavy  mass  concentration  in 
the  fuselage  of  modern  fighter  aircraft  results  in  the 
tendency  for  the  x-body  axis  to  become  perpendicular  to 
the  axis  of  rotation.  If  this  is  allowed  to  happen,  a 
transformation  will  occur  where  the  angular  velocity 
vector  is  changed  from  consisting  entirely  of  the  body 
axis  roll  rate  to  consisting  entirely  of  the  body  axis 
yaw  rate  which,  of  course,  is  a  flat  spin.  The  direction 
of  the  yaw  build  up  is  a  function  of  both  the  direction 
of  the  roll  and  whether  positive  or  negative  angle  of 
attack  flight  conditions  are  present.  Experience  has 
shown  that  the  transformation  will  not  occur  smoothly, 
and  at  some  point  when  the  build  up  of  yaw  rate  exceeds 
the  amount  that  can  be  offset  by  the  inherent  directional 
stability  of  the  aircraft,  the  aircraft  will  become 
unstable,  and  a  departure  from  controlled  flight  will 
result.  This  is  a  physical  explanation  for  departures 
from  controlled  flight  as  a  result  of  inertial  roll 
coupling. 
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The  amount  of  roll  coupling  shown  in  Figure  7.4  is  not 
large,  however,  and  no  serious  roll  coupling  problems 
should  be  experienced  at  those  particular  flight 
conditions. 

e.  Case  Humber  Five  (DAI=-29.2.  DEL=-10.2.  DROrQ. -30) 

The  initial  starting  point  in  Figures  7.6a  through 
7.5e  is  the  full  scale  aileron  deflection  position  from 
case  number  two.  A  cross  control  input  of  right  rudder 
against  full  left  aileron  deflection  with  stabilator 
deflection  held  at  -10.2  degrees  was  accomplished.  This 
was  an  attempt  to  further  aggravate  the  adverse  yaw 
problem  noted  in  Figure  7.2.  The  eigenvalues  of  the 
individual  solution  points  were  once  again  all  stable.  A 
slight  increase  in  angle  of  attack  is  noted  in  Figure 
7.5b,  and  this  can  be  attributed  to  the  roll  coupling 
phenomena  discussed  above.  Figures  7.5c  and  7.5d  show  the 
direct  relationship  between  yaw  rate  and  roll  rate,  and 
that  the  effects  of  a  full  scale  aileron  deflection 
against  the  rudder  movement  appear  to  have  little  effect 
bn  the  roll  authority  of  the  rudder  at  these  flight 
conditions.  The  sideslip  angle  and  yaw  rate  plots  on  the 
next  page  once  again  show  the  tendency  for  the  aircraft 
to  rotate  about  the  flight  path  vector  with  sideslip 
remaining  near  zero. 


In  an  attempt  to  further  investigate  the  roll 
coupling  phenomenon  noted  in  case  number  four,  an 
equilibrium  surface  for  rudder  deflection  angles  of  0  to 
-30  and  0  to  18  degrees  was  generated.  The  DRU=0 
position  in  Figures  7.6a  through  7.6e  corresponds  to  the 
full  scale  aileron  deflection  point  in  case  number  four. 
Unlike  the  positive  angle  of  attack  cross  control 
equilibrium  surface  discussed  in  case  number  five,  the 
steady  state  flight  trajectories  shown  in  Figure  7.6  are 
far  more  complicated.  This  is  due  primarily  to  the 
combined  effects  of  the  inertial  roll  coupling  already 
present  at  the  starting  flight  condition  and  the  somewhat 
counterintuitive  behavior  created  by  rudder  deflections 
when  flying  at  negative  angles  of  attack.  Before 
presenting  Figures  7.6a  through  7.6e,  the  reader  is  once 
again  reminded  that  the  aircraft  is  already  under  the 
effects  of  inertial  coupling,  and  therefore,  the  body 
axis  rotation  rates  (p,  q,  and  r)  should  be  considered  as 
components  of  angular  rotation  about  the  flight  path.  A 
consequence  of  neglecting  gravity  effects  in  this  model 
allows  the  flight  path  to  be  considered  level  at  all 
times  without  any  loss  of  generality.  During  moderate  to 
high  angle  of  attack  flight,  the  blockage  effects  created 
by  the  horizontal  tail  make  the  pitching  moment,  the 
rolling  moment  and  the  roll  damping  due  to  rudder  surface 
deflection  negligible. 
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However,  at  negative  angles  attack  these  effects  are  not 
negligible,  and  they  play  an  important  roll  in  the 
results  shown  in  Figure  7.6.  When  the  rudder  is 
deflected  during  negative  angle  of  attack  flight,  several 
interrelated  inputs  to  the  aircraft  dynamical  situation 
are  made.  Given  a  right  rudder  deflection  (DRU=0  to  DRD=- 
30),  the  following  inputs  to  the  system  occur: 

a.  right  yawing  moment 

b.  positive  pitching  moment  (independent  of  rudder 
deflection  direction) 

c.  left  rolling  moment  (considering  the  tail  as  a 
wing  and  the  rudder  as  an  aileron) 

d.  left  sideslip  angle 

e.  alteration  of  the  flow  fields  around  the  tail 
section  resulting  in  increased  roll  damping  during 
left  rolls. 

With  the  rudder  deflection  inputs  discussed  above  in 
mind,  refer  to  Figure  7.6.  Figure  7.6c  depicts  the  most 
counterintuitive  results  of  this  particular  run  and  will 
be  discussed  first.  At  a  negative  angle  of  attack 
(provided  that  is  also  negative),  a  right  rudder 

deflection  should  result  in  a  left  rolling  moment.  This 
is  due  to  the  rolling  moment  due  to  yaw  rate,  sideslip 
and  rudder  deflection  all  having  the  same  algebraic  sign 
as  the  rudder  deflection  when  flying  at  negative  angles 
of  attack.  This  clearly  is  not  the  situation  indicated 


in  Figure  7.6c. 
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The  roll  rate  magnitude  is  sharply  decreased  initially 
when  the  rudder  is  moved  from  0  to  -30  degrees.  This 
occurence  can  best  be  explained  by  the  change  in  the  air 
flow  about  the  tail  that  occurs  when  the  rudder  is 
deflected.  During  a  roll*  the  horisontal  tail  provides 
some  damping  effect  but  it  is  usually  minimised  by  the 
blockage  of  airflow  due  to  the  vertical  tail.  When  a 
right  rudder  deflection  is  made*  the  roll  damping  due  the 
stabilator  is  enhanced  for  rolls  to  the  left  only.  The 
mechanism  responsible  for  this  phenomenon  is  high  energy 
airflow  directed  by  the  rudders  onto  the  right  upper 
stabilator  surface.  This  dramatically  increases  the 
damping  capabilities  of  the  horisontal  tail  and  is 
responsible  for  the  oounter intuitive  behavior  seen  in 
Figure  7.6c.  During  rudder  roll  manuevers  at  positive 
angles  of  attack,  however*  this  is  not  a  problem  since  a 
right  rudder  deflection  will  increase  the  damping  of  left 
rolls  only  thereby  improving  the  capability  of  the  rudder 
to  check  roll  rates  while  at  the  same  time  not  decreasing 
the  roll  rate  generated  by  the  rudder  itself.  Figure 
7.6e  shows  that  the  sideslip  angle  changes  sign  at  a 
rudder  deflection  of  approximately  -8  degrees.  This 
corresponds  to  the  inflection  point  seen  in  the  roll  rate 


curve. 
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Ultimately,  the  negative  sideslip  angle  is  responsible 
for  producing  a  left  rolling  moment  strong  enough  to 
overcome  the  increased  roll  damping  due  to  rudder 
deflection.  When  considering  the  left  rudder  input  in 
this  situation,  the  damping  effect  will  not  be  present 
since  the  aircraft  is  rolling  to  the  left  in  this 
equilibrium  branch.  The  remaining  results  are  relatively 
straight  forward.  In  Figure  7.6a  the  right  rudder 
deflection  (DEL=0  to  DEL=-30)  clearly  shows  a  smooth 
increase  in  pitch  rate.  This  is  due  to  both  the  pitching 
moment  created  by  the  rudder  deflection  and  the  decreased 
roll  rate  seen  in  Figure  7.6c  which  causes  a  reduction  in 
the  negative  pitching  moment  due  to  inertial  coupling. 
During  a  left  rudder  deflection.  Figure  7.6c  shows  the 
roll  rate  to  eventually  decrease  as  the  rudder  deflection 
exceeds  10  degrees.  The  increasing  sideslip  angle  shown 
in  Figure  7.6e  is  the  driving  force  behind  that 
behavior.  In  Figure  7.6d  the  yaw  rate  decreases  when 
rudder  is  deflected  in  either  direction.  This  phenomenon 
is  more  easily  understood  by  breaking  the  yaw  rate  up 
into  two  separate  components,  one  component  is  due  to 
rudder  deflection  and  the  other  component  is  due  to 
inertial  coupling.  At  the  point  where  DRU=0,  the  yaw  rate 
shown  is  made  up  primarily  of  yaw  rate  due  to  inertial 
coupling.  A  left  rudder  deflection  naturally  tends  to 
deorease  the  yaw  rate. 


130 


In  the  case  of  the  right  rudder  deflection,  the  decrease 
in  yaw  rate  due  to  inertial  coupling  is  greater  than  the 
increase  in  yaw  rate  due  to  rudder  input  thus,  an  overall 
decrease  in  yaw  rate  is  demonstrated  initially.  The  two 
yaw  rate  components  appear  to  balance  each  other  out  as 
the  rudder  deflection  is  decreased  below  -20  degrees. 

The  changes  in  stability  noted  in  the  DRU=7  to 
DRU=18  degrees  range  shows  the  possibility  for  a  Hopf 
Bifurcation  to  a  limit  cycle  as  DRU  is  increase  past  7 
degrees  and  the  possibility  of  a  small  nonlinear  jump  at 
DRU=16  degrees.  The  numerical  convergence  was  excellent 
throughout  this  particular  run  including  at  the  point  in 
question.  This  fact  as  well  as  the  positions  of  the 
eigenvalue  in  question  before  and  after  the  unstable 
point  indicate  that  a  very  small  bifurcation  whose  limit 
points  were  stepped  over  by  the  computer  algorithm  may 
exist  in  that  area  of  the  curve. 

g.  Case  Number  Seven  <DAI=-29.2,  DRU=-29.2.  DEL-- 10. 2 

The  starting  point  for  this  equilibrium  surface 
comes  from  case  number  five  and  is  a  pitch  up  manuever 
with  full  cross  controls.  The  aileron  is  held  at  -29.2 
degrees  and  the  rudder  is  held  at  -29.2  corresponding  to 
full  loft  aileron  and  full  right  rudder. 


The  results  indicate  a  distinct  bifurcation  in  the 
equilibrium  surface  in  the  -15  to  -20  degrees  stabilator 
range.  In  this  case  the  figures  are  shown  one  per  page 
to  provide  a  better  view  of  the  results.  Some  step  size 
and  convergence  difficulties  were  encountered  in  the 
computation  of  the  portion  of  the  curve  connecting  the 
two  limit  points.  The  curves  were  generated  by  initially 
varying  the  stabilator  deflection  from  -10.2  to  -30. 
The  run  went  very  well  numerically  but  a  jump  in  both  the 
equilibrium  surface  and  the  locations  of  the  eigenvalues 
was  noted  at  approximately  -17  degrees  stabilator 
deflection.  The  run  was  repeated  with  a  smaller  step 
size  and  another  equilibrium  surface  was  generated  in  the 
other  direction  by  varying  stabilator  deflection  from  -20 
degrees  to  -10  degrees.  The  repeat  of  the  original  run 
continued  to  a  slightly  lower  stabilator  setting  before 
making  a  similar  jump  and  then  following  the  same 
trajectory  as  the  original  run.  The  reverse  direction 
run  followed  the  same  trajectory  as  the  original  run 
except  that  the  jump  was  made  at  a  significantly  higher 
stabilator  defleotion  than  the  point  where  a  jump  was 
made  during  the  first  run  thus  showing  a  definite 
hysteretic  pattern.  These  results  were  confirmed  by 
subsequent  runs,  all  of  which  showed  good  numerical 
convergence,  but  the  proper  step  size  to  continue  around 
the  limit  points  was  not  found. 
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This  indicates  that  some  problems  may  exist  with  the 
continuation  algorithm.  Although  an  actual  calculation  of 
the  eigenvalues  of  the  branch  connecting  the  two  limit 
points  was  not  computed,  the  behavior  of  the  trajectory 
along  that  branch  indicates  that  it  is  probably  an 
unstable  branch.  A  smooth  curve  was  fit  between  the  last 
confirmed  point  on  the  upper  and  lower  branches  to  give  a 
more  realistic  view  of  the  equilibrium  surface  in  that 
region.  The  hash  marks  just  proir  to  the  limit  points 
show  the  location  of  the  last  confirmed  solution  points 
prior  to  the  jump. 

Figure  7.7a  shows  the  equilibriuim  pitch  rate  versus 
stabilator  deflection  for  the  conditions  given  above.  The 
pitch  rate  shows  little  change  up  to  the  limit  point 
where  at  that  point  a  small  but  sharp  drop  in  pitch  rate 
is  predicted.  The  point  dropped  to  is  right  on  the  S/L 
stability  boundary  and  further  movement  of  the  stabilator 
places  the  aircraft  in  a  Hopf  Bifurcation  region  where 
limit  cycle  motion  may  be  experienced.  Given  the  type  of 
manuever  and  the  large  jump  in  roll  rate  that  occurs 
during  the  catastrophe  (shown  in  figure  7.7c),  a  wing 
rock  type  motion  is  a  solid  possibility.  Comparison  of 
Figure  7.7a  with  Figure  7.7b  shows  that  contrary  to  the 
other  fifth  order  surfaces  the  angle  of  attack  behavior 
at  the  bifurcation  point  is  opposite  of  the  pitch  rate 
behavior. 


A  small  but  sharp  increase  in  angle  of  attack  is 
predicted  at  DEL=-16.5.  Figures  7.7c  through  7.7e  all 
show  similar  nonlinear  jumps  with  the  largest  one 
occur ing  in  the  roll  rate  equilibrium  surface.  The  sharp 
jump  in  sideslip  angle  from  slightly  less  than  zero  to 
approximately  -4  degrees  corresponds  closely  to  the  large 
increase  in  roll  and  raw  rates  shown  in  Figures  7.7c  and 
7.7d.  In  each  of  the  plots  except  for  the  angle  of  attack 
equilibriuim  surface,  the  location  of  the  L  instability 
corresponds  to  an  area  where  the  state  variable  shown  is 
changing  sharply  with  stabilator  deflection.  Since  the 
catastropic  changes  shown  in  the  aircraft  states  must  be 
accompanied  by  a  similar  change  in  kinetic  energy,  I 
believe  that  a  limit  cycle  does  exist  in  the  region 
where  the  L  instability  occurs,  and  it  is  through  the 
limit  cycle  oscillations  that  the  kinetic  energy  is  built 
up  or  released  prior  to  or  after  the  nonlinear  jump. 
Whether  the  energy  is  increased  or  decreased  during  the 
limit  cycle  motion  is  dependent  on  the  direction  from 
which  the  bifurcation  region  is  approached.  As  cm 
example,  consider  the  case  of  increasing  stabilator 
deflection  from  DEL=-25  to  DEL=-10.  In  this  situation 
the  equilibrium  roll  rate  is  shown  to  begin  to  increase 
significantly  with  increasing  stabilator  deflection  at 
DEL=-20  which  is  also  the  point  where  the  L  instability 
is  encountered. 


As  the  aircraft  state  crosses  the  stability  boundary,  a 
limit  cycle  is  born.  The  limit  cycle  behavior  during  this 
phase  imparts  energy  to  the  system  and  roll  rate  is 
increase  until  the  eigenvalues  are  again  all  stable. 
This  corresponds  to  the  level  portion  of  the  upper  roll 
rate  branch  where  the  increase  in  roll  rate  stops.  The 
system  is  now  ripe  for  a  catastrophe  which  does  indeed 
occur  at  DEL=-15  The  energy  previously  gained  is  then 
released  from  the  system  during  the  catastrophe  so  that 
the  energy  level  is  once  again  at  the  level  prior  to 
entering  the  limit  cycle  phase.  The  aircraft  behavior 
during  the  catastrophe  cannot  be  predicted  by  the 
equilibrium  surface,  but  it  is  clear  that  it  will 
involve  the  dissapation  of  the  energy  built  up  as  the 
catastrophe  was  approached.  A  definite  hysteresis  effect 
can  be  seen  in  this  hypothesis  as  well  as  in  all  of  the 
equilibrium  surface  trajectories  for  this  case.  If  the 
bifurcation  is  approached  from  the  other  direction,  the 
energy  increase  occurs  suddenly  in  a  manner  very  similar 
to  a  departure  from  controlled  flight.  The  energy  gained 
during  the  catastrophic  jump  at  DEL=-16.5  is  dissipated 
gradually  by  the  limit  cycle  motion  as  the  stabilator 
deflection  angle  is  decreased  until  the  kinetic  energy  is 
returned  to  a  level  approximately  equal  to  the  level  of 
kinetic  energy  prior  to  the  catastrophe. 
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The  information  contained  in  these  equilibrium 

surface  plots  is  physically  substantiated  in  the  results 

reported  in  a  NASA  Technical  Memorandum  by  Grafton, 

et  al.  [8].  In  the  document  the  results  of  a  stability 

study  on  a  wind  tunnel  model  very  similar  to  the  F-15S 

are  reported.  In  that  report  the  canard  stalling  angle 

of  attack  is  shown  to  be  approximately  30  degrees.  This 

corresponds  to  the  angle  of  attack  where  the 

bifurcational  behavior  occurs.  At  that  angle  of  attack  a 

sharp  stabilizing  break  in  the  C  versus  angle  of  attack 

m 

curve  as  well  as  a  very  sharp  peak  in  the  versus  angle 

of  attack  curve  occurs.  This  explains  in  part  the 

behavior  of  the  angle  of  attack  and  pitch  rate  curves. 

The  change  in  sign  of  the  pitching  moment  stability 

derivative  at  that  point  can  account  for  the  sharp 

decrease  in  pitch  rate  that  is  seen  across  the 

catastrophe.  The  angle  of  attack  in  the  fifth  order 

model  is  a  function  of  both  pitch  rate  and  C  .  The  sharp 

z 

loss  of  lift  that  occurs  as  the  slope  of  the  lift  curve 
changes  signs  will  create  a  plunging  effect  along  the  z- 
body  axis  which  tends  to  increase  angle  of  attack  and  is 
most  likely  responsible  for  the  catastrophe  seen  in  the 
angle  of  attack  equilibrium  surface.  The  report  also 
states  that  at  an  angle  of  attack  of  30  degrees  the 
aircraft  "exhibits  highly  nonlinear  characteristics  with 
a  very  abrupt  unstable  variation  in  rolling  moment  over  a 
small  sideslip  range  near  0  sideslip  angle." 
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Theae  instabi 1 itiea  were  found  to  be  a  direct  result  of 
the  canard  configuration  and  disappeared  when  the  canard 
was  removed.  Additionally,  the  report  states  that  this 
configuration  exhibits  its  minimum  yaw  damping  at  an 
angle  of  attack  of  30  degrees  and  that  around  35  degrees 
angle  of  attack  there  is  "an  abrupt  decrease  in  roll 
damping  to  very  low  values."  This  is  precisely  the  type 
of  behavior  predicted  in  the  yaw  and  roll  rate 
equilibrium  surface  plots,  and  the  sideslip  plot  does 
confirm  that  the  aircraft  is  near  zero  sideslip  angle. 
And  finally,  the  report  also  indicates  that  this 
configuration  can  be  expected  to  be  susceptible  to  wing 
rock  motion  in  the  angle  of  attack  range  corresponding  to 
the  L  instability  region  which  adds  some  credence  to  the 
energy  transfer  hypothesis  stated  above.  These  findings 
correlate  very  closely  with  the  predictions  made  by  the 
bifurcation  and  catastrophe  analysis  shown  in  Figures 
7.7a  through  7.7e.  No  known  simulations  of  this  type  of 
configuration  in  this  or  any  other  similar  high  angle  of 
attack  flight  regime  have  been  accomplished.  Therefore, 
the  existence  of  the  limit  cycle  energy  transfer  model 
hypothesized  above  can  not  readily  be  established.  But 
there  is  clear  evidence  that  the  bifurcations!  behavior 
exhibited  in  Figures  7.7a  through  7.7b  does  exist,  and 
the  behavior  predicted  through  their  study  is  correct. 


As  was  discussed  above  only  two  eighth  order 
equilibrium  surfaces  were  calculated.  The  increased 
degrees  of  freedom  offered  by  the  eighth  order  model  not 
only  increased  the  complexity  of  the  results  but 
tremendously  increased  the  numerical  difficulty  of  the 
problem.  The  major  factor  contributing  to  the  bulk  of 
the  numerical  problems  appeared  to  be  the  extra  degree  of 
freedom  allowing  velocity  to  vary.  This  caused  the 
equilibrium  surface  calculations  with  thrust  to  be  very 
difficult  since  velocity  tended  towards  either  zero  or  a 
very  small  value  as  the  routine  searched  for  the  initial 
equilibriuim  points.  The  larger  the  thrust  the  faster 
velocity  would  tend  to  zero.  This  behavior  indicated 
conditions  typical  of  operating  an  aircraft  ‘behind  the 
power  curve* .  Therefore,  initial  conditions  of  level 
flight,  thrust  equal  to  drag  and  velocity  equal  to  a 
value  well  above  that  for  minimum  drag  were  selected. 
The  results  were  the  same  and  study  of  the  system  with 
thrust  was  discontinued.  This  presented  no  major 
problems  since  the  study  of  spin  type  motions  was  the 
main  goal  of  the  eighth  order  investigation. 
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Velocity  variations  did  continue  to  cause  some  problems, 
and  although  the  problems  were  not  as  severe  as  those 
present  with  thrust,  the  large  velocity  variations 
during  the  iterations  to  initial  solution  points  as  well 
as  new  points  during  the  continuation  process  ultimately 
resulted  in  only  (me  non-spin  equilibrium  surface  being 
computed.  The  one  that  was  found  could  not  be  continued 
to  a  stabilator  deflection  greater  than  -5  degrees.  This 
corresponded  to  an  angle  of  attack  of  approximately  29 
degrees.  Grafton  et  al  [8]  reports  that  the  canard  had 
serious  detrimental  effects  on  pitch  damping  at  angles  of 
attack  above  25  degrees  and  this  may  partially  explain 
the  problems  encountered  in  this  area.  Another  intersting 
aspect  is  that  while  the  fifth  order  system  was  mostly 
stable,  the  eighth  order  system  was  rarely  stable.  This 
will  be  discussed  more  thoroughly  in  the  next  section. 

a.  Case  Number  One  (DAI=Q.  DRU=0.  DEL=-5  to  -3Q) 

Case  one  was  a  pitch  up  manuever  with  aileron  and 
rudder  deflections  held  to  zero.  The  results  of  this  run 
are  presented  in  Figures  7.8a  through  7.8p.  Several 
plots  are  duplicated  for  comparison  purposes  and  those 
that  are  duplicated  are  given  the  same  title.  Some  plots 
are  repeated  in  an  amplified  scale  to  see  the  small  order 
variations  that  occur  in  the  equilibrium  surface. 
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Figures  7.8a  and  7.8b  show  the  pitch  angle  and  angle  of 
attack  equilibrium  surfaces.  These  plots  indicate  that 

the  aircraft  has  a  descent  angle  that  varies  from 
approximately  2  degrees  at  DEL=-5  to  a  descent  angle  of 
6  degrees  at  DEL=-30.  The  shape  of  the  two  curves  is 
nearly  identical  and  comparison  with  the  pitch  rate  plot 
in  Figure  7.8c  on  the  next  page  show  that  including  the 
gravity  coupling  terms  has  decoupled  the  pitch  rate 
behavior  from  the  angle  of  attack  behavior  and  at  the 
same  time  has  apparently  coupled  the  pitch  angle  behavior 
to  the  angle  of  attack.  A  zero  pitch  rate  is  realistic 
given  the  high  angle  of  attack  conditions  and  zero 
thrust.  Figure  7.8d  shows  an  amplified  version  of  the 
pitch  rate  equilibrium  surface  which  reveals  a  sharp 
.  change  in  behavior  as  DEL  is  decreased  below  -25  degrees. 
This  point  corresponds  to  an  angle  of  attack  of  40 
degrees  and  it  is  at  this  point  where  the  aerodynamic 
model  removes  the  drag  input  due  to  ailerons  by  setting 
K2  to  zero  (see  Appendix  D).  Additionally,  two  other 
impulsive  type  changes  are  made  to  the  system.  One  is 
riade  at  32  degrees  angle  of  attack  which  corresponds  to 
the  point  where  the  UUU  instability  transitions  to  the 
UUL  instabilty  and  the  other  is  at  an  angle  of  attack  of 
35  degrees  which  corresponds  to  the  point  where  the  UUL 
instability  transtions  to  an  L  instability. 
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Also,  the  first  impulsive  change,  discussed  above,  occurs 
at  the  point  where  transition  from  the  L  instability  to 
the  LL  instability  occurs..  Study  of  the  actual 
aerodynamic  data  shows  that  the  terms  involved  in  the 
impulsive  changes  were  decreased  to  zero  or  near  zero 
values  prior  to  the  angles  of  attack  given  above  to 
create  a  smooth  transition;  however,  it  is  hard  to 
believe  that  these  occurances  are  unrelated,  and  there 
may  be  a  problem  with  the  rate  of  change  of  those  values 
not  also  being  zero  at  the  point  where  they  are  included 
or  removed  from  the  aerodynamic  data.  Figure  7.8e  shows 
the  aircraft  velocity  equilibrium  surface,  and  the 
behavior  of  the  curve  is  what  should  be  expected  given 
the  flight  conditions.  Figures  7.8f  and  7.8g  show  the 
roll  and  yaw  rate  equilibrium  surfaces,  and  these  confirm 
the  symetrical  nature  of  the  aerodynamic  model.  Amplified 
versions  of  the  roll  and  yaw  surfaces  are  shown  on  the 
next  page,  and  they  reveal  the  same  type  of  behavior  seen 
in  the  previous  amplified  plot  as  well  as  a  definite 
coupling  between  the  two  states.  The  next  two  pages  show 
the  sideslip  angle  and  yaw  rate  plots  together  in  both 
the  non-amplified  and  the  amplified  versions.  The 
amplified  versions,  Figures  7.8k  and  7.81  show  nearly 
identical  behavior  between  yaw  rate  and  sideslip. 
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The  fact  that  the  sideslip  angle  and  yaw  rate  share  the 
same  algebraic  sign  is  an  indication  of  the  increased 
complexity  introduced  in  the  model  by  removing  the 
restriction  that  gravity  be  equal  to  zero.  In  this 
particular  case  a  very  slight  “dutch  roll”  type  behavior 
is  noted  which  is  a  type  of  motion  suppressed  in  the 
fifth  order  system.  Figures  7.8n,  f,  o  and  p  show  the 
relationship  between  roll  rate  and  bank  angle. 

b.  Case  Number  Two  (DEL=-21.  DRU=-24. 5,  DAI =-29. 2) 

This  equilibrium  surface  is  the  spin  surface 
discussed  earlier.  It  was  found  by  simply  trying  various 
estimated  spin  conditions  until  the  N ewt on -Raph s on  method 
was  able  to  converge  to  an  equlibrium  solution.  As  has 
been  mentioned  before,  no  spin  studies  have  been  done  on 
this  configuration.  The  trajectory  given  by  the 
equilibrium  surface  routine  is  quite  reasonable  and 
indicates  a  strong  possibility  for  recovery  at  full  high 
performance  aircraft  anti -spin  controls.  In  the  case  of 
this  particular  spin  which  is  a  left  spin,  the  anti-spin 
controls  are  full  deflection  of  left  aileron,  right 
rudder  and  negative  stabilator  where  the  aileron  is 
usually  considered  the  primary  spin  recovery  control 
surface.  The  type  of  spin  trajectory  indicated  is  a  flat 
spin  with  very  little  variation  in  each  of  the  state 
variables  except  yaw  rate  with  the  movement  of  aileron 
defleotion. 
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The  scales  used  for  the  Y  axis  in  the  plots  vary  from 
plot  to  plot  to  provide  the  ‘best*  view  of  the  individual 
surfaces.  Figures  7.9a  and  7.9b  show  the  relationship 
between  pitch  angle  and  angle  of  attack  as  aileron  is 
varied.  The  curves  are  nearly  identical  and  the  sum  of 
their  angles  as  measured  from  the  vertical  is  90  degrees 
indicating  a  vertical  descent  condition.  As  the  aileron 
deflection  nears  -30  degrees  both  of  the  curves  drop  off 
sharply  indicating  a  possibility  for  recovery.  Figure 
7.9c  shows  the  variation  of  yaw  rate  with  aileron  and  a 
very  large  increase  in  spin  rate  is  seen  as  the  left 
aileron  input  is  decreased.  Figure  7.9d  shows  the 
velocity  equilibrium  surface  and  that  is  compared  to  the 
angle  of  attack  surface.  The  velocity  shows  little 
variation  except  when  near  DAI =-30.  Figure  7.9e  shows 
the  sideslip  angle  equilibrium  surface  and  that  too  shows 
evidence  for  recovery  at  full  anti-spin  aileron  input. 
The  greatest  evidence  of  possible  spin  recovery  is 
indicated  in  Figure  7.9g  on  the  next  page.  In  this  figure 
the  pitch  rate  curve  actually  undergoes  a  large  change  in 
slope  as  the  anti-spin  aileron  input  is  increased. 
Figures  7.9i  and  7.9j  show  an  enlarged  view  of  the  pitch 
rate  and  velocity  surfaces  taken  from  an  equilibrium 
surface  run  accomplished  in  an  attempt  to  obtain  a 
surface  for  the  full  left  aileron  deflection. 
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Evidence  of  numerical  difficulties  is  shown  by  the  poor 
convergence  of  a  couple  of  the  points  near  the  end  of  the 
run. 

Ideally,  obtaining  a  surface  that  extends  all  the 
way  to  DAI =-30  degrees  as  well  as  increasing  the  other 
two  controls  to  full  anti-spin  deflections  would  help 
verify  the  proper  recovery  control  strategy,  but  serious 
numerical  difficulties  were  encountered  when  this  was 
attempted.  This  problem  stems  from  the  same  situation 
found  during  the  generation  of  the  fifth  order 
equilibrium  surface  in  the  neighborhood  of  limit  points 
for  case  number  seven.  Unfortunately,  in  the  case  of  the 
spin  surface  the  jump  that  probably  occurs  at 
approximately  DAI =-29. 5  degrees  is  quite  large  and  the 
algorithm  is  unable  to  converge  to  this  new  point. 
Overall,  the  equilibrium  spin  trajectory  depicted  in  this 
case  appears  to  be  correct  with  one  exception  that  will 
be  discussed  below.  The  plots  clearly  indicate  the 
possibility  of  a  limit  point  ocurring  between  DAI=-29  and 
DAI=-30  degrees.  If  that  is  the  case,  then  the  standard 
anti-spin  controls  should  create  a  jump  from  the  spin 
equilibrium  surface  shown  to  some  other  surface  that  is 
at  least  a  lower  energy  spin  surface  or  quite  possibly  a 
full  recovery  to  controlled  flight. 


A  noteworthy  point  indicating  that  this  surface  may 
not  predict  an  actual  stabilized  spin  condition  for  the 
aircraft  is  given  by  the  local  stability  analysis  being 
unstable  throughout.  This  indicates  that  the  surface  is 
an  unstable  equilibrium  branch,  and  the  likelihood  of  a 
similar  stable  equilibrium  branch  that  exists  nearby  is 
present.  A  few  solution  points  were  found  during  the 
study  indicating  that  a  lower  energy  stable  spin  surface 
may  indeed  exist,  but  numerical  problems  once  again 
prevented  the  complete  computation  of  the  branch 
associated  with  those  few  points. 

It  should  be  remembered  at  this  time  that  the  rotary 
balance  'data  used  may  not  be  as  accurate  as  is  necessary 
to  give  the  proper  stability  indications.  The  data  was 
good  enough  however  to  allow  the  calculation  of  a 
reasonable  spin  equilibrium  trajectory.  An  analysis  of 
spin  data  for  the  F-15  £45]  showed  that  the  F-15  should 
have  stabilized  spin  rates  in  the  150  to  170  deg/sec 
range,  an  angle  of  attaok  in  the  80  to  85  degree  range 
and  velocity  in  the  230  to  250  f/s  range.  These  values 
are  reasonably  close  to  those  given  by  the  spin 
equilibrim  surface  in  this  study.  The  addition  of  the 
canard  can  account  for  a  good  deal  of  the  difference  in 
the  stabilized  spin  modes  of  the  two  aircraft. 
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During  high  angle  of  attack  spin  conditions,  the  canard 
will  provide  a  destabilizing  pitch  up  increment  since  it 
is  above  its  stalling  angle  of  attack  and  therefore  is 
operating  in  the  region  of  negative  lift  curve  slope. 
This  fact  when  combined  with  the  fact  that  the  canard 
adds  additional  planform  area  forward  of  the  center  of 
gravity  indicates  that  the  F-15S  should  have  a  flatter 
stabilized  spin  mode  when  compared  to  the  original  F-15. 
An  aircraft  experiencing  a  relatively  flatter  spin  tends 
to  have  a  higher  yaw  rate  since  the  aircraft  center  of 
gravity  is  closer  to  the  spin  axis  and  the  other  rotation 
rates  are  decreased  to  near  zero  values.  Therefore,  as  a 
result  of  conservation  of  angular  momentum,  the  yaw  rate 
component  is  increased. 


The  numerical  simulation  portion  of  the  study  did 
not  yield  a  great  deal  of  useful  information.  The  model 
simulated  was  an  eighth  order  model  which  was  shown  in 
the  local  stability  analysis  to  be  unstable,  and  it  was 
indeed  shown  to  be  unstable  during  the  simulations  that 
will  be  presented  in  this  section. 
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The  instability  of  the  model  was  not  considered  to  be  a 
great  draw  back  initially  since  it  was  intended  to  use 
the  computer  aided  stability  system  during  the  bulk  of 
the  numerical  simulations.  Two  sets  of  two  different  but 
related  simulations  will  be  presented  here.  The  first 
set  is  at  a  relatively  high  angle  of  attack  without 
thrust,  and  the  second  set  is  at  a  lower  angle  of  attack. 
One  of  the  two  runs  in  the  second  set  is  with  thrust. 
The  scales  on  the  plots  are  not  constant. 


During  these  two  simulations,  the  initial  conditions 
were  identical  and  they  were  taken  from  the  conditions  in 
the  eighth  order  non-spin  equilibrium  branch  for  DEL=-8.7 
except  that  the  stabilator  was  set  at  -5.0  as  the  initial 
condition  for  the  simulation.  In  Figure  7.10  the  results 
of  a  smooth  decrease  in  stabilator  deflection  angle  is 
shown  and  in  Figure  7.11  a  smooth  increase  in  stabilator 
deflection  is  shown  up  to  where  DEL=+. 1  and  then  held 
constant.  Initially,  a  decrease  in  angle  of  attack  is 
expected  due  to  the  pitch  up  stabilator  setting  being 
less  than  that  required  for  equilibrium  conditions  and 
that  is  shown  to  be  the  case  in  both  of  the  runs. 


However,  from  that  point  stabilator  deflection  seems  to 
have  no  effect  on  the  aircraft.  The  control  inputs  were 
double  checked  and  appeared  to  be  correct.  The 
subroutine  that  generates  the  aerodynamic  coefficients  as 
a  function  of  the  control  surface  deflections  as  well  as 
a  number  of  other  variables  is  the  same  one  used 
successfully  in  the  fifth  and  eighth  order  equilibrium 
surface  calculations.  This  pointed  to  the  possibility 
that  the  aircraft  may  actually  have  very  limited 
stabilator  control  authority  in  the  situation  shown.  The 
next  obvious  step  was  to  determine  the  canard  behavior 
during  this  simulation.  The  canard  is  scheduled  as  a 
function  of  angle  of  attack  alone  and  the  schedule  is 
shown  below: 

DC  =  2.0  -  1.25  *  OC  ( 7 . 1 ) Th 
e  canard  deflection  range  is  -36  to  4-18  degrees  with 
positive  deflection  being  trailing  edge  down.  The  canard 
schedule  is  designed  to  unload  the  canard  at  high  angles 
of  attack  to  delay  the  adverse  effects  created  by  the 
canard  stall  and  thus  extend  the  usable  angle  of  attack 
range.  Therefore,  at  25  degrees  angle  of  attack  the 
canard  is  at  -29.25  degrees  deflection  and  increases  to  - 
10.5  at  an  angle  of  attack  of  10  degrees.  As  a  result  of 
this,  the  effects  of  the  canard  in  the  case  shown  in 
Figure  7.10a  tend  to  work  against  the  stabilator 
deflection. 
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The  nearly  identical  angle  of  attack  curves  indicate  that 
the  canard  may  have  more  control  authority  than  the 
stabilator  in  the  situation  shown  in  the  two  figures. 
Figures  7.10b  and  7.11b  show  identical,  near  zero 
behavior  for  roll  rate,  sideslip  angle  and  pitch  rate. 
This  verifies  the  symmetric  trajectory  behavior  predicted 
by  the  non-spin  eighth  order  equilibrium  surface.  Figure 
7.10c  shows  an  enlarged  plot  of  the  decrease  in  angle  of 
attack  plotted  against  stabilator  deflection. 

2.  Cases  Three  and  Four  ( DA 1=0.  DEL=-1.  DRU=0) 

Figures  7.12a  through  7.12d  show  the  results  of  a 
five  second  simulation  where  the  initial  state  was  taken 
from  the  DEL=-5  degrees  non-spin  equilibrium  branch.  The 
actual  stabilator  deflection  used  in  this  run  as  well  as 
case  four  was  DEL=-1 .  In  these  plots  a  large  decrease  in 
velocity  with  the  increase  in  angle  of  attack  is  shown  up 
to  the  point  where  the  angle  of  attack  reaches 
approximately  70  degrees.  At  this  point  it  is  apparent 
that  the  aircraft  is  developing  a  ‘sink  rate'  as  the 
velocity  is  again  increased  due  to  the  increasing 
component  in  the  z-body  axis  direction.  Figures  7.12c 
and  7.12d  show  a  slow  but  definite  instability  in  pitch 


rate  and  roll  rate. 
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The  sideslip  angle  (not  shown)  showed  a  smooth  increase 
from  0  up  to  +.77  degrees  at  56  degrees  angle  of  attack 
and  then  a  smooth  decrease  at  nearly  the  same  rate  to 
-.22  degrees  where  the  simulation  aborted  at  90  degrees 
angle  of  attack.  An  interesting  point  to  note  at  this 
time  is  that  55  degrees  is  the  angle  of  attack  where  the 
rotary  balance  data  is  included  in  the  aerodynamic  model, 
and  this  may  account  for  the  decrease  in  sideslip  angle 
seen  past  that  point.  The  final  simulation  was  similar  to 
case  three  except  an  initial  value  for  velocity  that 
should  have  placed  the  aircraft  or.  the  'front  side'  of 
the  power  curve  was  chosen.  The  initial  velocity  was  640 
f/s,  and  the  data  for  Mach=0.6  was  used.  Roll  rate,  yaw 
rate,  pitch  rate  and  sideslip  angle  were  all  set  to  zero 
and  remained  there  throughout  the  simulation. 
Additionally,  thrust  was  set  to  8500  lbs  in  this  case. 
Figures  7. 12a  and  7. 12b  show  the  velocity  and  angle  of 
attack  versus  time  curves.  The  curves  shown  in  these 
figures  are  very  similar  to  the  ones  in  Figures  7.12a  and 
7.12b  with  the  two  exceptions  that  a  longer  period  of 
time  was  required  in  case  four  to  reach  90  degrees  angle 
of  attack  and  the  velocity  range  traversed  in  case  four 
was  higher  than  case  three. 
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In  the  simulations  described  above,  the  eighth  order 
model  is  shown  to  be  definitely  unstable  as  is  predicted 
by  the  local  stability  analysis  in  eighth  order  case 
number  one.  The  lack  of  significant  lateral  motion 
predicted  by  the  equilibrium  surface  when  the  rudder  and 
ailerons  are  held  at  zero  degrees  deflection  is  also 
verified.  The  addition  of  thrust  appears  to  change  the 
time  scale  of  the  instability  but  not  the  qualitative 
nature  of  the  unstable  trajectory.  Little  more  than  this 
can  be  said  concerning  the  numerical  simulations 
themselves  or  their  comparison  to  the  equilibrium 
surfaces.  As  a  result  of  the  obvious  limits  to  the 
usefulness  of  numerical  simulations  without  the  computer 
aided  stability  system,  the  time  spent  on  the  numerical 
simulations  during  this  study  was  less  than  the  other 
areas  discussed  in  this  chapter. 

D.  Comparison  of  the  Fifth  and  Eighth  Order  Models 
In  this  section  the  significant  differences  and 
similarities  between  the  two  equilibrium  surface  models 
will  be  discussed.  The  main  areas  of  interest  in  this 
portion  of  the  discussion  involve  comparisons  of  the 
model  trajectories,  stability  and  numerical  behavior. 
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A  significant  difference  between  the  conditions  that 
the  two  investigations  were  accomplished  under  involves 
the  thrust  being  zero  in  the  case  of  the  eighth  order 
calculations  and  the  velocity  being  held  constant  in  the 
case  of  the  fifth  order  calculations.  Holding  velocity 
constant  gives  the  fifth  order  trajectory  a  certain 
amount  of  'implied  thrust'.  This  is  seen  in  the  angle  of 
attack  versus  stabilator  deflection  equilibrium  surface 
plots  of  the  two  models  shown  in  Figures  7.1b  and  7.8b. 
Both  of  the  models  show  a  maximum  angle  of  attack  near  40 
degrees  at  full  scale  stabilator  deflection,  but  the 
angle  of  attack  range  hat  the  fifth  order  equilibrium 
surface  shows  is  much  larger  than  the  eighth  order  plot. 
For  instance,  the  equilibrium  angle  of  attack  at  DEL=-5 
degrees  is  nearly  10  degrees  higher  in  the  eighth  order 
system  than  the  fifth  order  system.  This  is  a  direct 
effect  of  allowing  velocity  to  vary  in  the  eighth  order 
model.  As  discussed  above  the  pitch  rate  and  angle  of 
attack  are  decoupled  in  the  eighth  order  model  while  on 
the  other  hand,  the  pitch  angle  and  angle  of  attack 
appear  to  be  coupled  in  the  eighth  order  model.  In  both 
of  the  models  the  lack  of  lateral  motions  during  a  pitch 
up  manuever  is  shown,  and  this  was  confirmed  to  a  certain 
extent  in  the  numerical  simulations  discussed  above. 
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The  increased  degrees  of  freedom  allowed  in  the  eighth 
order  model  create  a  more  subtle  coupling  between  the 
lateral  terms.  This  is  seen  in  the  amplified  plots  shown 
in  Figure  7.8  while  that  type  of  behavior  was  not  present 
in  the  fifth  order  pitch  up  equilibrium  surface. 

In  terms  of  local  stability  behavior  the  differences 
are  quite  significant.  The  fifth  order  equilibrium 
points  tended  to  be  locally  stable  during  most  of  the 
investigation  while  just  the  opposite  was  noted  in  the 
eighth  order  equilibrium  points.  The  numerical 
simulations  that  were  accomplished  combined  with 
discussions  with  personnel  involved  with  the  F-15S 
project  at  Wright-Patterson  AFB,  Ohio  support  the  eighth 
order  system  as  having  the  correct  stability 
information.  The  behavior  of  both  the  fifth  order  model 
and  the  equilibrium  curves  indicate  that  the  fifth  order 
stability  calculations  are  correct  for  the  model  used 
(additionally,  the  same  routine  was  used  to  calculate  the 
eigenvalues  for  both  of  the  equilibrium  surface 
routines),  but  it  is  apparent  that  the  released  degrees 
of  freedom  in  the  eighth  order  model  contain  the 
critical  destabilizing  components  found  in  the  actual 
aircraft.  Consequently,  it  is  believed  that  the  fifth 
order  model  is  inadequate  to  make  stability  predictions 


on  the  F-15S. 
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On  the  other  hand,  however,  the  equilibrium  trajectories 
revealed  by  the  fifth  order  model  seems  to  be  at  least 
qualitatively  correct,  and  the  bifurcational  behavior 
found  in  case  number  seven  does  explain  the  existence  of 
stability  problems  revealed  during  earlier  research  on 
the  configuration.  These  seemingly  contradictory  results 
can  be  reconciled  by  refering  to  the  theoretical 
discussion  on  the  concepts  of  bifurcations  and 

catastrophes  in  Chapter  4.  Applying  the  concept  proposed 
by  the  center  manifold  theorem  to  the  local  stability 
analysis  of  the  eighth  order  system  shows  that  there  is  a 
group  of  three  to  five  eigenvalues  that  remain  in  the 
neighborhood  of  the  imaginary  axis  with  the  other 
eigenvalues  far  removed  either  to  the  right  or  left  of 
the  imaginary  axis.  Comparison  of  these  eigenvalues  with 
the  ones  present  in  the  fifth  order  case  show  those  three 
to  five  eigenvalues  also  to  be  present  in  the  fifth  order 
system.  Therefore,  it  appears  that  the  reduced  order 
system  as  determined  by  the  center  manifold  theorem  is 
contained  within  the  fifth  order  system  used  in  this 
investigation,  and  thus,  the  complete  bifurcational 
behavior  of  the  aircraft  should  be  revealed  in  the 
behavior  of  the  fifth  order  model.  The  added  degrees  of 
freedom  in  the  eighth  order  model  do  tend  to  make  the 
aircraft  unstable,  but  this  behavior  is  linear  in  nature 
and  is  therefore  irrelevant  to  the  bifurcational  behavior 
of  the  system. 
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The  complete  behavior  of  the  system,  according  to  Arnol'd 
[37]  can  be  found  from  the  fifth  order  system  by 
'suspending'  the  linear  unstable  and  stable  parts  present 
only  in  the  eighth  order  system  to  the  nonlinear  fifth 
order  system.  The  resulting  system  will  be  identical  to 
the  eighth  order  system  in  this  study  except  the 
bifurcational  behivior  will  have  been  found  by  studying 
the  much  simpler  fifth  order  system.  This  accounts  for 
the  differences  in  local  stability  behavior  of  the  two 
systems,  and  at  the  same  time  shows  the  close 
relationship  between  the  two  models  from  a  bifurcational 
behavior  point  of  view.  It  is  clear  now  that,  contrary  to 
initial  indications,  the  fifth  order  system  can  indeed 
reveal  accurate  information  regarding  the  behavior  of  the 
complete  aircraft  model,  and  provides  a  means  to  at  least 
partially  sidestep  the  numerical  problems  encounted  in 
the  eighth  order  system  by  studying  the  bifurcational 
behavior  of  the  fifth  order  system  and  then  generalizing 
the  results  to  the  eighth  order  system.  This  would  reduce 
the  number  of  eighth  order  equilibrium  surface  runs  which 
as  will  be  shown  in  the  next  section  could  provide  a 
significant  improvement  in  the  efficiency  and  success  of 
similar  studies. 

Numerically  speaking,  the  differences  between  the 
two  models  is  very  large. 
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After  fine  tuning  the  algorithm  and  including  a  number  of 
computer  time  saving  steps,  some  of  which  are  discussed 
in  chapter  six,  the  fifth  order  equlibrium  surface 
calculation  routine  was  working  quite  well.  By  the  end 
of  the  study,  none  of  the  partial  differentiation  or 
convergence  problems  found  in  the  eighth  order  routine 
were  present  in  the  fifth  order  routine,  and  the  computer 
time  per  solution  point  was  reduced  to  under  30  seconds 
during  some  of  the  fifth  order  runs.  The  eighth  order 
system,  on  the  other  hand,  was  not  nearly  so  easy  to  work 
with.  The  difficulties  were  isolated  for  the  most  part 
to  the  large  variation  caused  in  each  of  the  states  by 
small  changes  in  velocity.  In  an  attempt  to  solve  this 
problem  weighting  factors  that  weighted  the  relative 
importance  of  having  the  time  derivatives  of  each  of  the 
state  variables  be  zero  at  the  solution  point  were  used. 
Both  a  relatively  small  value  and  a  relatively  large 
value  for  the  velocity  time  derivative  weighting  factor 
were  attempted  with  no  real  success.  Problems  were  even 
encountered  during  attempts  to  converge  to  a  previously 
established  equilibrium  point.  Finally,  both  of  the 
routines  appeared  to  have  difficulties  continuing  the 
equilibrium  surface  around  limit  points.  This  problem 
may  be  a  result  of  either  improper  continuation  step  size 
or  improper  weighting  of  the  columns  in  the  augmented 
system  Jacobian  matrix  (see  Chapter  5); 


however,  a  number  of  runs  were  attempted  where  these 
factors  were  varied  and  still  little  success  was  found. 
This  tends  to  indicate  that  the  method  used  may  need  to 
be  improved  to  be  useful  on  a  more  complex  or  unstable 
aerodynamic  model. 
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CHAPTER  VIII 

SUMMARY  OF  SIGNIFICANT  FINDINGS 

In  this  chapter  the  more  significant  results 
discovered  during  the  course  of  this  study  are 
summarized.  The  order  in  which  they  are  presented  here 
is  the  same  relative  order  in  which  they  were  originally 
presented  in  the  above  discussion  to  aid  the  reader  in 
locating  the  more  indepth  discussions  on  these  findings 
located  in  Chapter  7. 

During  the  fifth  order  equilibrium  surface 
investigation,  three  of  the  surfaces  yielded  behavior 
that  predict  notable  problem  areas  for  the  F-15S.  The 
first  of  these  findings  is  that  the  aircraft  is 
susceptible  to  severe  adverse  yaw  problems  in  the  25 
degree  angle  of  attack  range,  and  this  results  in  rolling 
moments  being  produced  that  are  opposite  of  the  aileron 
control  deflection  direction.  This  phenomenon  is 
demonstrated  in  Figures  7.2a  and  7. 2d  where  a  left 
aileron  input  is  shown  to  generate  a  smooth  yawing  motion 
to  the  right  which  in  turn  results  in  a  smooth  rolling 
motion  to  the  right  at  a  stabilator  deflection  of  -10.2 
degrees. 

The  second  area  of  interest  stems  from  the  roll 
coupling  behavior  seen  at  an  angle  of  attack  of 
approximately  -5  degrees. 
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The  inertial  roll  coupling  forces  on  the  aircraft  when 
combined  with  additional  control  inputs  create  highly 
nonlinear,  counter-intuitive  behavior  in  the  equilibrium 
trajectories.  The  first  evidence  of  roll  coupling 
behavior  is  shown  in  Figure  7.4  where  pitch  rate  and 
angle  of  attack  are  seen  to  decrease  with  increasing  roll 
rate  to  the  left  while  the  yaw  rate  smoothly  increases. 
The  unusual  behavior  described  above  occurs  when  a  rudder 
control  movement  is  made  after  reaching  an  aileron 
deflection  of  -29.2  degrees  (full  left  position)  in 
Figure  7.4.  The  -29.2  degrees  aileron  position  in  figure 
7.4  and  the  0  degrees  rudder  position  in  figure  7.6 
correspond  to  the  same  point  on  an  aileron  versus  rudder 
equilibrium  surface.  The  highly  nonlinear, 
counterintuitive  behavior  mentioned  above  occurs  as  a 
result  of  an  interaction  between  the  inertial  roll 
coupling  forces  present  prior  to  the  rudder  input,  the 
forces  that  are  a  direct  result  of  the  control 
deflections  and  the  forces  created  as  a  result  of  a 
variation  in  the  flow  fields  about  the  aircraft.  In  the 
case  shown  in  Figure  7.6  a  rudder  deflection  to  the  right 
is  seen  to  decrease  the  left  roll  rate  by  causing  a 
variation  of  the  flow  field  about  the  aircraft  tail 
section  which  in  turn  causes  a  large  increase  in  roll 
damping.  The  lower  roll  rate  decreases  the  effects  of 
roll  coupling  which  when  combined  with  the  other  forces 
create  the  nonlinearities  seen  as  well  as  aircraft  state 
changes  that  are  opposite  of  the  expected  behavior. 
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These  interactions  and  the  mechanisms  causing  them  are 
more  completely  discussed  in  Chapter  7.  The  discovery  of 
this  unusual  behavior  does  underscore  the  usefulness  of 
the  equilibrium  surfaces  in  studying  aircraft  behavior 
during  unusual  control  sequences.  By  studying  aircraft 
trajectories  from  an  equilibriuim  surface,  the  complex 
transient  behavior  that  would  normally  accompany  these 
flight  conditions  is  removed  so  that  the  underlying 
causes  of  the  behavior  can  be  studied  in  a  more  direct 
manner. 

The  third  and  most  important  finding  during  the 
fifth  order  equilibrium  surface  investigation  is  the 
discovery  of  an  area  with  dirtinct  bifurcational 
behavior.  Figures  7.7a  through  7.7e  show  the  evidence 
of  catastrophic  behavior  in  the  aircraft  states  during 
stabilator  movements  in  the  -20  to  -15  degrees  stabilator 
deflection  range  while  holding  full  left  aileron  and  full 
right  rudder.  The  nonlinear  jump  behavior  predicted  by 
this  finding  is  substantiated  in  a  NASA  Technical 
Memorandum  by  Grafton  et  al.  [8]  in  which  the  stability 
of  an  aircraft  very  similar  to  the  F-15S  is  studied 
through  wind  tunnel  tests.  In  the  NASA  document  the 
existence  of  abrupt  changes  in  the  aircraft  stability 
derivatives  at  the  angle  of  attack  where  the 
bifurcational  behavior  is  seen  to  occur  in  Figure  7.7  is 
reported  but  not  explained.  The  occurance  of  the 
behavior  reported  in  the  NASA  document  is  accurately 
predicted  by  catastrophe  theory. 
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Additionally,  a  mechanism  based  on  a  Hopf  Bifurcation  to 
limit  cycle  motion  is  proposed  to  explain  the  physical 
transfer  of  kinetic  energy  that  must  occur  during  the 
nonlinear  jumps.  The  basis  for  this  mechanism  is 
partially  supported  by  one  of  the  predictions  made  in  the 
NASA  Technical  Memorandum.  According  to  that  document, 
the  origin  of  the  nonlinear  behavior  is  traced  to  the 
addition  of  the  canard  to  the  original  F-15 
configuration. 

During  the  eighth  order  equilibrium  surface 
investigation  a  spin  trajectory  equilibrium  surface  was 
discovered  for  variations  in  aileron  deflection.  The 
values  of  the  state  variables  correspond  reasonably  close 
to  known  stabilised  spin  modes  for  the  original  F-15 
configuration.  The  spin  mode  found  for  the  F-15S  is 
flatter  and  has  a  higher  rotation  rate  than  that  typical 
of  the  F-15,  but  these  differences  in  spin  behavior  can 
be  attributed  to  the  addition  of  the  canard  to  the  basic 
F-15  configuration.  Additionally,  the  spin  equilibrium 
surface  indicates  a  significant  possibility  for  a  jump 
from  the  spin  surface  shown  to  either  a  lower  energy  spin 
or  possibly  to  non-spin  flight  conditions  by  using  a  full 
anti-spin  aileron  deflection.  One  significant  detraction 
from  the  reliability  of  this  spin  mode  prediction  is  the 
locally  unstable  eigenvalues  present  at  each  of  the 
equilibrium  points  indicating  that  the  spin  surface  shown 
is  a  series  of  unstable  equilibrium  points.  For  a  more 
detailed  discussion  on  this  point  see  Chapter  7. 
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Due  to  the  unstable  nature  of  the  aircraft  without 
an  operational  computer  aided  stability  system,  very 
little  useful  information  was  gained  from  the  numerical 
simulations.  One  significant  point  that  was  noted, 
however,  was  that  during  a  power  off  stall,  stabilator 
deflection  appeared  to  have  little  bearing  on  the 
aircraft  behavior.  This  effect  is  shown  in  Figures  7.10 
and  7.11.  The  canard  behavior  was  determined  to  have  an 
input  that  opposes  pitch  up  stabilator  deflections  at  the 
conditions  simulated,  and  that  may  explain  the  behavior 
noted. 

The  comparisons  between  the  eighth  and  fifth  order 
local  stability  results  show  a  large  disparity.  The 
evidence  indicates  that  the  fifth  order  model  does  not 
accurately  predict  the  local  stability  behavior  of  the 
actual  aircraft,  but  the  trajectory  and  bifurcation 
information  gained  during  the  fifth  order  equilibrium 
surface  investigation  was  shown  to  be  useful  in 
predicting  aircraft  behavior.  The  discussion  on  this 
topic  explains  that  the  large  difference  in  the  local 
stability  of  the  two  models  is  not  necessarily 
significant  in  terms  of  studying  the  bifurcational 
behavior  of  the  aircraft  since  the  additional  three  roots 
in  the  eighth  order  model  tended  to  be  far  removed  either 
to  the  left  or  to  the  right  of  the  imaginary  axis. 
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Refering  to  the  center  manifold  theorem  (Chapter  4) 
showed  that  the  reduced  order  system  for  the  F-15S  is, 
for  the  most  part,  contained  within  the  fifth  order 
system,  and  therefore,  the  bifurcational  complete 
behavior  of  the  aircraft  could  be  studied  by  studying  the 
fifth  order  system.  The  complete  aircraft  behavior  is 
then  recovered  by  generalizing  the  fifth  order  nonlinear 
behavior  to  the  eighth  order  system  in  the  manner 
discussed  by  Arnol'd  [37].  Once  this  is  done,  the 
complete  behavior  of  the  aircraft  including  the  linear 
instabilities  seen  in  the  eighth  order  model  as  well  as 
the  nonlinear  bifurcational  behavior  found  in  the  fifth 
order  model  will  be  seen,  and  of  course,  the  end  result 
will  be  locally  unstable  wherever  either  model  predicts 
local  instabilities. 

In  some  cases  the  findings  reported  above  are 
primarily  of  a  qualatative  nature,  however,  each  of 
these  as  well  as  those  findings  based  on  more 
quantitative  results  are  shown  to  be  valid  based  on  the 
documentation  available  and  the  evidence  presented  in  the 
results  of  this  study.  Further  investigation  of  the 
areas  listed  above  through  the  use  of  more  specific 
techniques  such  as  those  discussed  in  Chapter  3  is 
warranted,  and  additionally.  the  application  of  this 
technique  to  this  area  of  flight  mechanics  research 
should  be  continued  so  that  a  better  understanding  of 
this  and  other  high  angle  of  attack  flight  phenomena  can 
be  gained. 
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CHAPTER  IX 

CONCLUSIONS  AND  RECOMMEND ATA IONS 

A,  Summary 

During  this  research,  a  method  for  studying  aircraft 
flight  mechanics  through  the  use  of  bifurcation  and 
catastrophe  theories  was  studied  and  implemented  into  the 
computer  facilities  available  at  the  Mas sachus setts 
Institute  of  Technology.  The  study  was  based  primarily 
on  the  generation  of  equilibrium  surfaces  for  a  fifth  and 
an  eighth  order  aircraft  model.  The  aircraft 
configuration  studied  was  an  F-15  with  a  canard  mounted 
on  the  engine  inlet  shelves.  The  aircraft  is  known  as 
the  F-15  STOL  Demonstrator  and  is  currently  being 
developed  and  studied  by  McDonnell  Douglas  and  the  Air 
Force  Wright  Aeronautical  Laboratory  at  Wright-Patterson 
AFB.  The  aircraft  is  inherently  unstable  and  requires  a 
computer  aided  stability  system  for  normal  operation.  The 
results  in  this  study  are  for  the  aircraft  without  the 
computer  aided  stability  system.  Seven  fifth  order  and 
two  eighth  order  equilibrium  surfaces  were  successfully 
computed,  and  those  results  are  all  reported  in  this 
document.  Numerical  simulations  were  attempted,  but  due 
to  the  unstable  nature  of  the  aircraft  little  information 


was  gained  from  them. 
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A  summary  of  the  significant  findings  of  this  research  is 
contained  in  Chapter  8,  and  a  more  indepth  discussion  of 
all  of  the  results  as  well  as  those  significant  findings 
discussed  in  Chapter  8  are  found  in  Chapter  7. 

Overall,  the  research  was  a  successful  effort,  a 
good  deal  of  the  objectives  set  were  to  at  least  some 
degree  realized.  One  major  shortcoming  of  the  research 
was  the  inability  to  determine  more  tangible  results 
based  on  the  bifurcational  behavior  of  the  aircraft.  Time 
and  .logistical  limitations  were  directly  responsible  for 
the  second  major  limitation  which  was  the  failure  to  make 
any  studies  of  the  aircraft  with  its  computer  aided 
stability  system  intact.  The  amount  of  work  required  to 
make  the  software  package  compatible  with  the  computer 
facilities  available  and  to  implement  it  into  the 
software  already  developed  was  much  greater  than 
anticipated.  At  this  point,  I  would  like  to  thank  Julie 
Wolf  for  the  time  and  effort  she  expended  working  on  the 
implementation  of  the  flight  control  system.  It  was  not 
by  any  means  an  easy  task. 
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B.  Recommendations 

There  are  a  number  of  areas  in  which  this  research 
has  led  to  the  identification  of  suggestions  for  possible 
improvements  and  further  research.  A  few  of  these  areas 
are  aircraft  design  improvements,  improvements  in  the 
numerical  methods  used  and  areas  where  further  study  may 
be  beneficial.  These  points,  of  course,  are  listed  next. 

...  Suggestions for  General  Aircraft  Modifications 

The  adverse  yaw  problem  noted  in  the  fifth  order 
case  number  two  is  not  an  uncommon  occurence,  and  I  am 
sure  that  this  aspect  of  the  aircraft  behavior  has  been 
taken  into  account  in  the  aircraft’s  flight  control 
system.  If  this  problem  has  not  been  taken  into  account, 
then  as  a  minimum,  an  aileron  rudder  interconnect  system 
that  is  scheduled  as  a  function  of  angle  of  attack  should 
be  implemented.  The  lack  of  stabilator  authority  during 
the  power  off  post  stall  situation  discussed  in  the 
numerical  simulation  section  of  Chapter  7  indicates  that 
scheduling  the  symmetrical  deflection  of  the  canard  as  a 
function  of  angle  of  attack  and  Mach  number  alone  may  be 
a  bit  restrictive  in  post  stall  situations.  Finally, 
since  the  possibility  for  spin  recovery  shown  in  Figure 
7.9  does  not  occur  until  near  the  full  scale  aileron 
deflection,  extra  aileron  control  deflection  capability 
may  be  a  possible  means  for  insuring  the  control 
authority  to  recover  from  spin  situations. 
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This  extra  capability  could  be  included  whenever  the 
computer  aided  stability  system  senses  spin  type 
conditions . 

^Suggestions  for -Improvement  _of_ the  Numerical  Method 

The  efficiency  of  the  continuation  routines  was  a 
major  problem  during  a  good  deal  of  the  research. 
Improvements  were  made  in  the  algorithms  and  those  are 
discussed  in  the  text,  but  overall  the  method  is  still 
somewhat  time  consuming.  The  key  to  faster  operation  is 
to  decrease  the  number  of  derivatives  that  must  be 
calculated.  This  means  reducing  the  Jacobian 
calculations  or  ideally  eliminating  them  altogether.  The 
technique  of  only  calculating  the  Jacobian  when  the 
aircraft  state  travels  a  specified  distance  form  the  last 
time  the  Jacobian  was  calculated  worked  very  well  for  the 
fifth  order  system  and  made  some  substantial  improvements 
in  the  run  time  of  the  routine.  No  real  improvement  was 
noted  in  the  efficiency  of  the  eighth  order  routine. 
Another  helpful  technique  used  was  to  ensure  that  when 
the  aerodynamic  coefficients  are  computed  during  the 
Jacobian  calculation  phase.  only  the  terms  that  are 
varied  during  that  specific  call  to  the  aerodynamic 
coefficient  routine  are  recalculated.  This  too  was 
helpful  in  decreasing  run  time. 
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These  points,  however,  are  merely  fine  tuning  the 
efficiency  of  a  method  whose  cost  will  eventually  become 
prohibitive  as  the  order  and  complexity  of  the  system  is 
increased.  Ultimately,  a  derivative  free  algorithm  would 
give  the  greatest  margin  of  improvement  in  this  area. 

The  numerical  problems  noted  in  the  neighborhood  of 
the  limit  points  during  the  study  indicates  that  some 
improvement  is  needed  in  either  the  implementation  of  the 
continuation  method  used  or  the  actual  method  itself. 

C. -Suggestions  for  Further  Research 

The  most  obvious  research  step  to  take  next  is  to 
continue  the  research  on  the  F-15S  with  the  computer 
aided  stability  system  operational.  This  would 
undoubtedly  improve  the  usefulness  of  the  numerical 
simulations,  but  more  importantly,  the  model  being 
studied  will  then  closely  resemble  the  actual  aircraft 
being  developed  and  thus  the  information  revealed  can  be 
directly  used  to  further  the  development  of  the 
aircraft.  Another  very  interesting  area  to  study  will  be 
the  role  that  the  computer  aided  stability  system  will 
play  in  the  generation  of  the  equilibrium  surfaces.  It 
is  not  clear  whether  or  not  the  flight  control  system 
will  make  changes  in  either  the  local  stability  behavior 
of  the  solution  points  or  the  bifurcational  behavior  of 
the  aircraft. 


On  the  other  hand,  it  can  be  said  with  some  confidence 
that  since  the  flight  control  system  is  a  linear  full 
state  feedback  type  control  system  it  should  not  be  able 
to  effectively  deal  with  most  of  the  types  of  nonlinear 
instabilities  that  occur  as  a  result  of  bifurcational 
behavior.  The  comparison  of  a  computer  simulated 
trajectory  where  the  flight  control  positions  are  updated 
40  times  a  second  and  a  similar  equilibrium  surface  where 
the  flight  control  system  is  only  updated  on  each 
iteration  to  a  new  solution  point  would  be  most 
interesting.  The  freedom  that  the  flight  control  system 
has  to  vary  the  controls  between  iterations  may  make 
convergence  to  an  equilbrium  solution  point  very 
difficult,  or  it  may  actually  accelerate  the  convergence 
since  it  too  is  seeking  an  equilibrium  for  the  cockpit 
control  setting  given  to  it.  The  soil  for  research  in 
this  area  is  quite  fertile  and  is  virtually  unbroken. 

The  Hopf  Bifurcation  energy  transfer  model 
hypothesized  in  the  discussion  of  case  number  seven  has 
the  potential  to  warrant  further  study  by  continued 
equilibrium  surface  calculations  and  numerical 
simulations.  The  end  result  of  this  type  of  research 
would  be  the  development  of  ways  to  monitor  energy  build 
ups  of  the  type  discussed  in  case  number  seven  so  that 
impending  catastrophes  could  be  predicted  and  eventually 
prevented . 
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Refer ing  to  the  center  manifold  theorem  (Chapter  4) 
showed  that  the  reduced  order  system  for  the  F-15S  is, 
for  the  most  part,  contained  within  the  fifth  order 
system,  and  therefore,  the  bifurcational  complete 
behavior  of  the  aircraft  could  be  studied  by  studying  the 
fifth  order  system.  The  complete  aircraft  behavior  is 
then  recovered  by  generalizing  the  fifth  order  nonlinear 
behavior  to  the  eighth  order  system  in  the  manner 
discussed  by  Arnol'd  [37].  Once  this  is  done,  the 
complete  behavior  of  the  aircraft  including  the  linear 
instabilities  seen  in  the  eighth  order  model  as  well  as 
the  nonlinear  bifurcational  behavior  found  in  the  fifth 
order  model  will  be  seen,  and  of  course,  the  end  result 
will  be  locally  unstable  wherever  either  model  predicts 
local  instabilities. 

In  some  cases  the  findings  reported  above  are 
primarily  of  a  qualatative  nature,  however,  each  of 
these  as  well  as  those  findings  based*  on  more 
quantitative  results  are  shown  to  be  valid  based  on  the 
documentation  available  and  the  evidence  presented  in  the 
results  of  this  study.  Further  investigation  of  the 
areas  listed  above  through  the  use  of  more  specific 
techniques  such  as  those  discussed  in  Chapter  3  is 
warranted,  and  additionally,  the  application  of  this 
technique  to  this  area  of  flight  mechanics  research 
should  be  continued  so  that  a  better  understanding  of 
this  and  other  high  angle  of  attack  flight  phenomena  can 
be  gained. 
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APPENDIX  A 

F-15S  PHYSICAL  CHARACTERISTICS 
As  was  explained  in  chapter  two,  the  F-15S  is,  in 
nearly  all  aspects,  physically  identical  to  the  F-15B. 
It  is  an  F-15B  (serial  number  71-0290)  modified  to 
include  the  canards  and  the  two-dimensional  thrust 
vectoring  nozzles.  A  comprehensive  diagram  of  the 
planview  of  the  aircraft  is  shown  in  figure  A. 1.  The 
following  table  lists  the  precise  physical  charact¬ 
eristics  of  the  aircraft. 

TABLE  A. 1 


Takeoff  gross  weight  (clean 
internal  fuel) 

Ixx 

iyy 

Izz 

1X2 


configuration,  5718  lbs  of 
37,794  lbs  9 
25,938  slug-ft  , 
185,287  slug-ft“ 
206,359  slug— f 
-2,543  slug-ftz 


Ming 


Area  (reference) 

Area  (actual) 

Span 

Aspect  Ratio 
Taper  Ratio 
Sweep  (leading  edge) 
Dihedral 


Root  (BL  0)  NACA 

BL  77.0  NACA 

BL  155.0  NACA 

BL  224.73  NACA 

Tip  NACA 

Incidence 
Twist 

Modified  Conical 
Aileron  Area 
Aileron  Travel 


64A006.6 

64Ax05.9 

64Ax04.6 

64A203.5 

64A203.0 


Camber 


Flaperon  Area  (not  used) 
Flaperon  Travel  (not  used) 


608.00  sq  ft 
599.39  sq  ft 
42.808  ft 
3.01 
0.25 

45  degrees 
-1  degree 

301.5  in  actual 
226.0  in  actual 

149.6  in  actual 
94.0  in  actual 
68.3  in  actual 

None 

None 

CL/D=0.  3 

26.48  sq  ft 

±30  degrees 

35.84  sq  ft 

0  to  30  deg.  down 
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Stabilator 

Area  (reference) 

Area  (actual) 

Aspect  Ratio 

Taper  Ratio 

Sweep  ( leading  edge ) 

Dihedral 
Airfoil  &  -Chord 

Root  (BL  0)  NACA  0005.5-64 
BL  90.0  NACA  0003.5-64 

Tip  NACA  0002.5-64 

Stabilator  Travel 


120.00  sq  ft  total 
111.36  sq  ft  total 
2.05 
0.34 

50  degrees 
0 

137.2  in  reference 
117.9  in  reference 
46.5  in  reference 
4-20  to  -30  degrees 


Vertical  Tails 

Effective  Area 
Aspect  Ratio 
Taper  Ratio 
Sweep  (leading  edge) 

Airfoil  &  Chord 

Root  NACA  0005.0-64 

Tip  NACA  0003.5-64 

Rudder  Area 
Rudder  Travel 


125.2  sq  ft  total 

1.7 

0.27 

47 . 16  degrees 

115.0  in  actual 
30.6  in  actual 
19.94  sq  ft  total 
±30  degrees 


Canard 

Area 

Aspect  Ratio 
Taper  Ratio 
Span  (exposed) 

Sweep  (leading  edge) 

Dihedral 

Airfoil.  A-Chord 

Root  NACA  65A  w/sharp  LE  6% 
Tip  NACA  65A  w/sharp  LE  2X 
Canard  Travel  (Trailing  edge) 


88.1  sq  ft 

2.44 

0.46 

14.67  ft 

47. 16  degrees 

+20  degrees 

98.70  in 
45.45  in 

30  up,  20  dn  (deg) 


Since  nozzle  effects  and  configurations  other  than 
clean  were  not  considered  in  this  study,  details 
pertaining  to  those  areas  have  been  omitted  from  table 
A. 1  and  from  figure  A. 1  (shown  on  the  next  page). 


Pig.  A. 1  F-15S  Planform 
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APPENDIX  B 
ROTARY  BALANCE  DATA 

B. 1  General  Remarks 

Standard  wind  tunnel  aerodynamic  coefficient  data  is 
normally  composed  of  two  different  types  of  data,  static 
data  and  forced  oscillation  data.  The  fact  that  these 
types  of  data  neglect  the  complex  rotational  flows  that 
develop  around  an  aircraft  in  spin  type  motion  was  not  a 
factor  in  the  earlier  days  of  aircraft  development. 
Those  earlier  combat  aircraft,  in  the  same  manner  as  most 
general  aviation  aircraft  today,  could  rapidly  enter  and 
recover  from  a  spin  with  use  of  standard  aerosurface 
controls.  This  fact  when  combined  with  the  fact  that  the 
rotary  balance,  (RB),  data  for  those  earlier  aircraft 
showed  little  variation  with  changes  in  omega,  the  RB 
parameter  (to  be  discussed  in  detail  later),  allowed  the 
earlier  fighter  aircraft  to  be  successfully  modeled 
throughout  their  flight  envelope  without  using  rotary 
balance  data  (29).  This  is  not  the  case  with  modern 
fighter  aircraft.  Due  to  their  slender  bodies  and  low 
aspect  ratio  wings,  modern  fighter  aircraft  are  not 
easily  recovered  from  spin  type  situations,  and  the 
variation  of  their  RB  data  is  much  higher  than  their 
predecessors. 
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Although  currently  there  is  still  much  debate  over  the 
proper  implementation  of  RB  data  into  the  aerodynamic 
model,  it  is  clear  that  high  aoa  spin  type  conditions  can 
not  be  properly  modeled  without  RB  data. 

B. 2  Implementation  of  the  RB  Data 
Rotary  Balance  wind  tunnel  tests  are  conducted  by 
rotating  the  wind  tunnel  model  about  the  relative  wind 
with  the  axis  of  rotation  being  through  the  CG.  There 
are  two  major  points  of  contention  in  the  RB  data 
debate.  The  first  is  when  to  include  the  RB  data  in  the 
aerodynamic  model  and  the  second  is  how  the  RB  parameter, 
omega,  should  be  modelled.  A  consensus  has  been  found, 
more  or  less  on  the  first  point  by  the  physical 
limitations  of  the  current  testing  apparatus  and 
procedures.  Due  to  the  nature  of  the  wind  tunnel  test 
setup  for  rotary  balance  testing,  data  is  not  usually 
taken  at  angles  of  attack  below  55  degrees.  This  limit 
arises  from  the  fact  that  the  sting  is  connected  to  the 
center  of  gravity  through  the  top  of  the  aircraft. 
Therefore,  rotating  the  aircraft  at  lower  angles  of 
attack  becomes  impossible  due  to  geometric  and 
aerodynamic  interference  between  the  tail  and  the  sting. 


Sideslip  angles  greater  than  +10  degrees  are  also  not 
normally  attempted  due  to  the  large  hole  that  must  be 
made  in  the  top  of  the  aircraft  to  allow  for  that 
particular  degree  of  freedom  and  still  have  the  sting 
connected  to  the  CG.  The  sise  of  the  hole  necesary  for 
testing  at  sideslip  angles  greater  than  ±10  would  distort 
the  flow  fields  to  the  point  of  destroying  the  accuracy 
of  the  results.  A  combination  of  these  circumstances  and 
observations  of  actual  aircraft  spin  behavior  make  an  aoa 
of  55  degrees  the  usual  RB  data  input  aoa.  The  RB  data  is 
usually  retained  in  the  aerodynamic  coefficient  models 
until  the  aoa  drops  below  30  degrees.  The  removal  of  the 
RB  data  at  an  aoa  of  30  degrees  is  somewhat  arbitrary, 
but  it  is  based  on  the  fact  that  the  aoa  threshold  for 
spin  entry  is  normally  higher  than  the  aoa  at  which 
rotational  motion  still  occurs  during  a  recovery  from  a 
spin.  Additionally,  some  spin  test  setups  allow  the 
gathering  of  data  at  angles  of  attack  as  low  as  30 
degrees  by  off-setting  the  aircraft  from  the  axis  of 
rotation.  If  this  type  of  data  is  unavailable,  then  the 
data  for  aoa=55  degrees  is  used  when  the  aoa  drops  below 
55  degrees. 

The  second  problem  area,  selecting  the  correct 
definition  for  omega,  is  a  much  more  difficult  task.  It 
is  possible  to  choose  several  definitions  that  have  merit 
at  one  particular  condition  but  fail  at  others. 
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2  2  2  2 

For  instance,  omega  =  p  +  q  +  r  or  even  omega=r  would 
adequately  describe  the  stabilized  flat  spin  condition 
since  p,q<<r  during  a  flat  spin,  but  neither  definition 
is  very  adequate  in  oscillatory  type  spin  conditions. 
This  is  a  major  drawback  of  most  of  the  models  considered 
and  arises  from  the  faot  that  the  RB  data  is  composed  of 
much  of  the  same  components  as  the  forced  oscillation 
data.  Therefore,  considerable  redundancy  is  introduced 
into  the  aerodynamic  coefficients  by  this  method.  This 
problem  is  normally  dealt  with  by  using  adjusted  values 
for  p,  q,  and  r  when  the  RB  data  is  in  use.  The  RB  data 
also  has  a  static  component  that  is  evident  by  the  non¬ 
zero  value  of  the  RB  coefficients  when  omega=0.  This  is 
removed  by  simply  subtracting  that  value  from  each  of  the 
RB  coefficients  calculated  so  that  RB=0  when  omega=0. 


M  •V1  i-i  A  »■»  nTWco  dRfCTQyK 


Since  no  RB  data  on  the  F-15S  has  been  collected  at 
the  present  time,  a  set  of  RB  data  for  the  F-15S  was 
approximated  by  adjusting  F-15  RB  data  to  account  for 
the  addition  of  the  canard.  The  F-15  RB  data  used  was 
presented  in  a  series  of  component  build-up  plots  that 
represented  the  various  changes  the  addition  of  specific 
components  had  on  the  rotary  balance  parameter, 
(omega*b}/(2*V)  (non-dimensional  form)  [44,45]. 
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This  type  of  data  presentation  made  it  possible  to 
estimate  the  effects  that  an  addition  of  the  canard  would 
have  on  the  RB  data  by  comparison  of  the  F-15  aerodynamic 
RB  data  for  omega=0  to  the  F-1SS  static  aerodynamic  data 
and  using  some  simple  calculations  to  find  the  relative 
effect  the  addition  of  the  canard  would  have  on  the  data 
based  on  the  effect  created  by  the  addition  of  the  other 
aerodynamic  surfaces.  Due  to  the  complex  nature  of  the 
flows,  no  estimations  were  made  for  data  other  than  zero 
sideslip.  A  discussion  of  the  process  used  to  develop 
the  F-15S  RB  data  is  given  below. 

The  Rotary  Balance  data  can  be  thought  of  as  a 
combination  of  static  and  forced  oscillation  data  with 
the  exception  that  rotational  flow  fields  are  allowed  to 
develop.  As  a  first  cut  approximation  the  values  for  the 
changes  in  static  forces  are  added  in  as  a  constant 
bias.  These  are  found  by  comparison  of  the  omega=0  RB 
data  for  the  F-15  to  the  corresponding  static  coefficient 
values  for  the  F-15S.  The  difference  is  the  static  bias. 
This  is  a  relatively  unimportant  point  for  our 
application,  however,  since  the  static  values  are  removed 
from  the  RB  data  to  avoid  duplication  of  the  static  force 
coefficients. 
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In  each  of  the  following  cases  where  adjustments  beyond 
the  changes  in  the  static  forces  are  required,  the 
adjustment  is  found  by  using  the  points  discussed  below 
to  determine  the  relative  effectiveness  of  the  canard  as 
compared  to  the  appropriate  tail  section  surface.  The  RB 
data  is  then  adjusted  accordingly  based  on  that 
particular  tail  section's  build-up  plot.  These  are  plots 
that  show  the  effect  that  the  addition  of  that  particular 
control  surface  alone  has  on  the  RB  data.  In  the  cases 
involving  the  longitudnal  coefficients  (Cm,  Cx  and  Cz) 
few  additional  adjustments  were  required.  Since  the 
canard's  stalling  aoa  is  equal  to  approximately  32 
degrees,  the  corrections  involved  only  simple 
calculations  of  the  canard  volume  versus  the  horizontal 
tail  volume  in  the  case  of  the  Cm  and  Cz,  and  a 
comparison  of  the  exposed  frontal  surface  area  of  the 
canard  versus  the  horizontal  tail  in  the  case  of  Cx 
(taking  into  account  the  blockage  and  downwash  effects  of 
the  wing  on  the  horizontal  tail).  In  the  calculation  of  * 
Cn  the  20  degree  dihedral  angle  of  the  canard  made  it 
necessary  to  make  some  adjustments  beyond  the  simple 
static  force  changes  discussed  above.  These  amounted  to 
recognition  of  the  destabilizing  contribution  of  the 
increased  vertical  surface  area  forward  of  the  CG  due  to 
the  canard,  and  a  comparison  of  the  yawing  moment  forces 
created  by  the  canard  with  those  that  would  be  created  by 
the  vertical  tails  without  the  stabilator. 
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The  canard  tends  to  be  destabilizing  in  yaw  at  zero  to 
low  rotation  rates  but  has  a  stabilizing  influence  at 
higher  rotation  rates  while  the  vertical  tail  is 
stabilizing  at  both  zero  and  high  rotation  rates.  This 
points  towards  the  probability  of  easier  spin  entries  but 
lower  stabilized  spin  rates  for  the  F-15S  when  compared 
to  the  original  F-15.  In  both  the  yawing  moment  and  the 
sideforce  coefficient  calculations  consideration  was 
given  to  the  effects  of  blockage  due  to  the  fuselage  on 
the  canard's  ability  to  influence  of  those  coefficients. 
Finally,  although  the  canard  has  some  significant 
influences  on  Cl,  it  loses  its  effectiveness  in  this  area 
when  placed  beyond  its  stalling  angle  of  attack,  and 
therefore,  only  small  adjustments  to  Cl  are  necessary. 
These  are  due  to  the  increased  stability  and  roll  damping 
offered  by  the  dihedral  angle  and  the  horizontal  surface 
area  of  the  canard. 

Although  the  above  adjustments  may  appear  crude  on 
the  surface,  it  should  be  remembered  that  any  data  is 
better  than  no  data  at  all  when  attempting  to  model  high 
aoa  spin  type  motions.  This  data  is  a  first 
approximation  to  the  actual  RB  data  and  is  therefore 
qualitatively  accurate  in  that  the  shift  of  the  data  is 
in  the  correct  direction,  but  the  amount  of  the  shift  may 
not  be  correct. 
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B.4  Application  of  the  RB  Data  to  the  F-1SS 
As  discussed  above  the  static  components  have  been 
removed  from  the  data  that  the  RB  coefficient  is  zero 
when  omega— 0 .  The  data  has  also  been  adjusted  so  that 

the  first  derivative  of  each  RB  coefficient  with  respect 
to  alpha  and  omega  is  0  when  aoa=55  and  30  degrees.  This 
prevents  an  impulsive  input  to  the  system  at  the  points 
where  the  data  is  to  be  included  and  removed  from  the 
system  of  aerodynamic  coefficient  equations.  Since 
sideslip  effects  are  not  in  the  data*  a  rotational 
parameter  was  selected  that  described  the  aircraft 
rotation  as  a  projection  of  the  total  aircraft  angular 
velocity  vector  onto  the  XZ  plane.  Thus, 

omega  =  p*cos(aoa)  4-  r*sin(aoa),  (B.4. 1) 

and  in  order  to  prevent  redundancies  between  the  forced 
oscillation  terms  and  the  RB  data  the  following 


substitutions  are  made  while  the  RB  data  is  being  used: 

P=P*sin(aoa)  (B.4. 2) 

r=r*cos(aoa)  (B.4. 3) 
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APPENDIX  C 

POINCARE-LIAPUNOV  STABILITY  THEOREM 

The  Poincare -Liapunov  Stability  Theorem  is  a  very 
simple  but  powerful  theorem.  It  has  been  the  basis  for 
the  successful  linearisation  of  nonlinear  systems  for 
many  years,  and  it  is  also  gives  major  support  to  the 
nonlinear  analysis  theorems  discussed  in  chapter  4.  In 
fact,  the  Center  Manifold  Theorem  is  merely  an 
extension  of  the  Pioncare-Liapunav  Theorem.  The 
Pioncare-Liapunov  Theorem  is  stated  below  without 
proof. 


"The  local  stability  of  a  nonlinear  system  at  a 
particular  operating  point  can  always  be  characterised 
in  the  same  way  as  the  equivalent  linear  system  that  is 
valid  in  the  neighborhood  of  that  operating  point 
except  in  those  cases  where  the  linearised  system  has 
eigenvalues  with  zero  real  parts." 
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APPENDIX  D 

F-15S  AIRCRAFT  MODELS 

This  appendix  contains  the  two  aircraft  models  used 
in  the  BACTM  software  and  the  model  used  in  the  numerical 
simulations.  The  BACTM  models  consist  of  a  fifth  and  an 
eighth  order  model  used  in  equilibrium  surface  calculations. 
A  separate  eighth  order  model  is  used  in  the  numerical 
simulation  phase. 

D.l  Aerodynamic  Force  Representation 
All  of  the  equations  of  motion  use  six  main  aerodynamic 
coefficients  that  are  calculated  in  the  exact  same  manner. 
This  is  true  for  cases  when  the  computer  aided  stability 
system  is  operational  and  when  it  is  not  operational.  The 
F-15S  data  is  represented  in  terms  of  CL,  CD,  Cy,  Cl,  Cm 
and  Cn.  The  lift  and  drag  coefficients  are  in  stability 
axis  coordinates  while  the  remaining  coefficients  are  in 
fixed  body  axis  coordinates.  The  lift  and  drag  coefficients 
are  transformed  into  body  axis  coordinates,  Cx  and  Cz,  for 
use  in  the  equations  of  motion.  For  clarity,  the  actual 
fortran  code  for  those  equations  is  used  in  the  descriptions 
shown  here. 

CL  ■  CLO  ♦  .5* (CLCANL+CLCANR)  ♦  CLFLX  ♦ 

.5*FLX1* (CLSTBL+CLSTBR)  ♦  DCLBET  (D.1.1) 

CD  *  CDO  ♦  .5* (CDCANL+CDGANR)  ♦  .5* (CDSTBL+CDSTBR)  ♦ 

K2 *DCDAI L * AVDA  ♦  CDRUD  +  DCDPCF  (D.l. 2) 
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CY1  =  CYO‘FLXNY*SBE  +  (1 . -FLXNY) *FLXCYO*BE  + 

DCYAIL*FLXAI *DA  +  CYDSTB*DHD*FLX6  ♦ 

.5* (CYRUDL*DRL+CYRUDR*DRR) *FLX5  +  CYCAN‘SBE  ♦ 
DCYDC*DCD  ♦  CYFLX  +  K7*DCYASY  (D.1.3) 

CLL1=  CLLO*FLXRM*SBE  ♦  DCLLP*P*B2V  ♦  DCLLR*R*B2V  ♦ 

DCLLAI *FLXAI *DA  ♦  CLLDST*FLX2*DHD  ♦ 

.5* (CLLRDL*DRL  +  CLLRDR*DRR) *FLX4  ♦  CLLCAN*SBE 
♦  DCLLDC*DCD  ♦  CLLFLX  (D.1.4) 

CM1=CM0  ♦  .5* (CMCANL+CMCANR)  «■  CMFLX  +  .5*FLX1* (CMSTBL 
+CMSTBR)  +  DCMAIL*AVDA  +  CMRUD  ♦  RAMP3*DCMDS2 
+  DCMQ*Q*C2V  ♦  .5* (CMSDCL+CMSDCR) *FLX1  «• 

AVDHD*CMCF  (D.1.5) 

CN1=  CNO*FLXNY*SBE  ♦  (1.0-FLXNY)  *FLXYM*BE  ♦  DCNP*P*B2V  ♦ 
DDCNP*K6*P*B2V  +  DCNR*R*B2V  ♦  DCNAIL*FLXAI *DA  + 
CNDSTB*FLX2  *DHD  +  .5* (CNRUDR ‘DRR+CNRUDL  *DRL) »FLX3 
+  CNCAN*SBE  «•  DCNDC'DCD  ♦  CNFLX  ♦  DCNASY*K7  (D.1.6) 

The  terms  used  In  the  above  set  of  equations  are  defined  as 

follows: 

P,  Q  and  R  are  roll,  pitch  and  yaw  rates  respectively, 

BE  is  the  sideslip  angle,  beta, 

B2V  and  C2V  are  span/(2*VT)  and  chord/ (2*VT)  respectively 

VT  stands  for  true  airspeed 

DA  is  the  aileron  surface  deflection. 

DCL  and  DCR  are  left  and  right  canard  surface  deflections. 

DHL  and  DHR  are  left  and  right  stabilator  surface  deflections. 

DRL  and  DRR  are  left  and  right  rudder  surface  deflections. 

DCD  and  DHD  are  differential  canard  and  rudder  surface 

deflections  and  are  defined  as: 

DCD=DCL-DCR  and 
DHD=DHR-DHL. 

DC,  DH  and  DR  are  the  symmetrical  canard,  stabilator  and 

rudder  deflections,  respectively  and  are  defined  as: 
DC= (DCL  ♦  DCR)/2 
DH=  (DHL  ♦  DHR)/2 
DR= (DRL  ♦  DRR)/2. 

A  prefix  of  AV  found  in  terms  such  as  AVDA  indicates  absolute 
value  of  the  term  following  the  AV. 

SBE  is  +1  whenever  beta  is  0  or  positive  and  -1  when  beta 
is  negative. 


The  individual  components  from  the  above  coefficient 
equations  are  grouped  together  below  into  terms  as  they 
appear  in  the  equations  and  the  overall  effect  of  each  term 
is  defined. 
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LIFT 

CLO 

.5* (CLCANL+CLCANR) 

CLFLX 

.5‘FLX1* (CLSTBL+CLSTBR) 

DCLBET 

DRAG 

CDO 

.5* (CDCANL+CDCANR) 

.5* (CDSTBL+CDSTBR) 

K2*DCDAIL‘AVDA 

CDRUD 

DCDPCF 


TERMS 

Basic  lift 
Lift  due  to  canard 
Canard  flexibility  effects 
Lift  due  to  stabilator 
Lift  due  to  beta 

TERMS 

Basic  drag 
Drag  due  to  canards 
Drag  due  to  stabilator 
Drag  due  to  ailerons 
Drag  due  to  rudder 
Per  formance  corr .  factor 


SIDE 

CYO*FLXNY*SB£ 

(l.-FLXNY) *FLXCYO*BE 

DCYAIL*FLXAI *DA 

CYDSTB*DHD*FLX6 

.5* (CYRUDL*DRL+CYRUDR*DRR) *1 

CYCAN*SBE 

DCYDC*DCD 

CYFLX 

K7*DCYASY 


:e  terms 

Side  force  due  to  beta 
Flexibility  effects 
Side  force  due  to  ailerons 
Side  force  due  to  stab. 

Side  force  due  to  rudder 
Side  force  due  to  sym.  can. 
Side  force  due  to  dif.  can. 
Canard  flexibility  effects 
High  alpha  asymmetry  effects 


CLLO*FLXRM*SBE 
DCLLP*P*B2V 
DCLLR*R*B2V 
DCLLAI *FLXAI *DA 
CLLDST*FLX2*DHD 


ROLLING  MOMENT  TERMS 

Rolling  mom. 
Rolling  mom. 
Rolling  mom. 
Rolling  mom. 
Rolling  mom. 


.5* (CLLRDL*DRL+CLLRDR*DRR) *FLX4  Rolling  mom. 
CLLCAN* SBE  Rolling  mom . 
DCLLDC*DCD  Rolling  mom. 


due  to  beta 
due  to  P 
due  to  R 
due  to  aileron 
due  to  stab, 
due  to  rudder 
due  to  sym.  can. 
due  to  dif.  can. 


CLLFLX 


Canard  flexibility  effects 


PITCHING 

CMO 

.5* (CMCANL+CMCANR) 

CMFLX 

.5*FLX1* (CMSTBL+CMSTBR) 

DCMAIL*AVDA 

CMRUD 

RAMP3*DCMDS2 

DCMQ*Q*C2V 

. 5* (CMSDCL+CMSDCR) *FLX1 
AVDHD*CMCF 


MOMENT  TERMS 

Basic  pitching  mom. 

Pitching  mom.  due  to  canard 
Canard  flexibility  effects 
Pitching  mom.  due  to  stab. 
Pitching  mom.  due  to  ailerons 
Pitching  mom.  due  to  rudder 
High  alpha  effects 
Pitching  mom.  due  to  Q 
Effect  of  canards  on  stab. 
Pitching  mom.  due  to  dif.  stab. 
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YAWING  MOMENT  TERMS 
CNO*FLXNY*SBE  Yawing  mom. 

(1.0-FLXNY) *FLXYM*BE  Flexibility 

(DCNP*P+DDCNP*K6)*P*B2V  Yawing  mom. 

DCNR*R*B2V  Yawing  mom. 

DCNAIL*FLXAI  *DA  Yawing  mom. 

CNDSTB * FLX2 *DHD  Yawino  mom. 

.5* (CNRUDR*DRR+CNRUDL*DRL) *FLX3  Yawing  mom. 

CNCAN*SBE  Yawing  mom. 

DCNDC*DCD  Yawing  mom. 

CNFLX 
DCNASY*K7 


due  to  beta 
effects 
due  to  P 
due  to 
due  to 
due  to 
due  to 
due  to 
due  to 


R 

aileron 
stab, 
rudder 
sym.  can. 
dlf.  can. 


Canard  flexibility  effects 
High  alpha  asymmetry  effects 


The  following  equations  show  the  transformations 


made  to  convert  the  above  stabillty/body  axis  set  of 
equations  into  a  body  axis  only  set  of  equations. 

CX=  CL*Sin (Alpha)  -  CD*Cos (alpha)  ♦  IROTCP*DCXRB  (D.1.7) 

CY=  CY1  ♦  IROTGP *DCYRB  (D.1.8) 

CZ-  -CD*Sln (alpha)  -  CL* Cos (Alpha)  ♦  IROTGP *DCZRB  (D.1.9) 

CLL=-CY1* (WLCG-WLREFCL) / (Span* 12 . )  +  CLL1  + 

IROTGP  *DCLLRB  (D.1.10) 

CM=  CX* (WLCG-WLREFCM) / (Chord* 12.)  -  CZ* (CG-CMCGR)  + 

CM1  +  IROTGP *DCMRB  (D.1.11) 

CN=  CY1*  (CG-CNCGR)  *CBAR/BB  ♦  CN1  +  IROTGP *DCNRB  (D.1.12) 


In  the  above  equations.  (D.1.7) -(D.1.12) ,  CG  is  the  Center 
of  Gravity  (measured  horizontally  in  %  of  chord) .  WLCG  is 
the  Water  Line  Center  of  Gravity  (measured  vertically  in 
Inches).  The  values,  WLREFCM,  WLREFCL,  CMCGR  AND  CNCGR, 
describe  the  location  of  the  origin  of  the  original  axis 
system  which  causes  moments  about  the  CG  in  the  new  body 
axis  system.  The  values,  DCXRB,  DCYRB.  DCZRB,  DCLLRB, 
DCMRB  AND  DCNRB.  are  the  rotary  balance  data  additions  for 
the  indicated  coefficients. 
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The  term,  IROTGP,  is  a  switching  parameter  that  is  0  whenever 
alpha  is  below  30  degrees.  It  is  switched  to  1  when  alpha 
increases  above  55  degrees  and  then  is  returned  to  0  when 
alpha  decreases  below  30  degrees.  See  appendix  B  for 
complete  details  on  the  rotary  balance  data. 

Finally,  one  last  adjustment  is  made  before  the 
above  coefficients  are  ready  to  be  used  in  the  equations 
of  motion.  If  thrust  is  to  be  Included,  then  CX  is  adjusted 
as  shown  below 

CX  s  CX  ♦  Thrust/ (.5*rho*VT** 2 *S)  (D.1.13) 

where  rho  is  the  atmospheric  density  and  S  is  the  wing 
planform  surface  area.  Figure  D.l  shows  the  sign  conven¬ 
tion  used  in  the  aerodynamic  model.  Note  that  the  aileron 
deflection  sign  convention  is  non-standard. 


Fiq.  D.l  F-15S  Sign  Conventions 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  OF  STANDARDS-1963-A 
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E_2  Equations  of  Motion 

As  discussed  above,  there  are  three  major  versions 
of  the  equations  of  motion  used  in  this  study.  The 
first  of  these  to  be  discussed  here  is  the  fifth  order 
set  used  for  equilibrium  surface  calculations.  This 
set  is  most  accurate  in  the  non-spin  flight  regime,  and 
it  is  derived  by  setting  gravity  to  zero  and  holding 
velocity  constant.  These  simplifications  reduce  the 
order  of  the  system  by  three.  Letting  gravity  be  sero 
decouples  the  bank  angle  and  pitch  angle  variables  from 
the  remaining  equations  of  motion,  and  therefore,  those 
terms  do  not  contribute  to  the  dynamics  of  the  system. 
It  has  been  shown  that  neglecting  gravity  does  not 
seriously  effect  the  modelling  of  the  flight  dynamics 
in  the  low  aoa  flight  regime,  but  the  assumption  of  g=0 
does  deteriorate  the  accuracy  when  modelling  high  aoa 
spin  type  behavior  [28,29].  Letting  velocity  be  con¬ 
stant  implies  that  the  thrust  available  is  always  equal 
to  the  thrust  required  which  is,  of  course,  unrealistic 
'  at  high  angles  of  attack.  The  restriction  of  velocity 
to  a  constant  value  also  has  the  effect  of  suppressing 
oscillations  in  the  other  varibles.  This  is  confirmed 
by  the  fact  that  wing  rock  and  other  oscillatory  type 
motions  are  not  readily  found  in  simulations  where 
velocity  is  not  allowed  to  vary  [28,29]. 
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Below  is  the  fifth  order  set  of  equations  used  in  this 
study  [28]. 


(D.2.1) 


<0.2. 2) 


♦ 


i&c 

mV 


cose 


(0.2.3) 


(0.2.4) 


■[ 

the  other  equilibrium  surface  set  of  equations 
also  uses  the  wind  axis  variables,  alpha  and  beta,  but 
it  also  includes  the  velocity,  pitch  angle  and  roll 
angle  equations  that  were  removed  to  oreate  the  fifth 
order  set  of  equations. 


sine 


-  p  si 


(0.2.5) 


na 
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It  is  one  of  several  sets  of  equations  that  can  fully 
describe  six  degrees  of  freedom  aircraft  motion.  Using 
the  wind  axis  variables  does  result  in  relatively 
complicated  equations  for  alpha  and  beta,  but  these 
variables  are  neoessary  since  they  contain  the  most 
meaningful  information  when  studying  high  aoa  spin  type 
motion.  By  calculating  them  directly  the  time  and 
effort  of  repeatedly  recalculating  their  values  as 
auxilliary  equations  is  saved.  The  eighth  order 
equilibrium  surface  set  of  equations  is  shown  below 
[6,28]. 


(D. 2.6) 


(D.2.7) 


(D. 2.8) 


(D.2.9) 
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-[(!■ c,  -  ?  si"°) s(ne  *  r]cos“ 
4(slCy  +  V  COS9Sin*)COS8 

*[(mV  cz  4  ?  cosecos^sine  -  pjs 


(D.2.10) 


-  f  sine)cosacos6  ♦  (&  cy  ♦  $  cosesln#) 
C2  ♦  $  cosocos^sinacosB 


sins 

(D.2.11) 


e  «  q  cos#  -  r  sin# 

#  *  P  ♦  Q  tanesin#  +  r  tanecos# 


(D.2.12) 

(D.2*13) 


The  set  of  equations  of  notion  selected  for  the 
numerical  simulations  is  also  an  eighth  order  set,  but 
in  this  case  the  fixed  body  axis  variables  (uf  v  and  w) 
were  chosen  over  their  wind  axis  counterparts  (alpha, 
beta  and  velocity).  This  sel -action  was  motivated  by 
the  clean  and  simple  form  of  the  body  axis  set  of 
equations.  It  was  felt  that  these  equations  would  make 
the  simulation  package  more  readily  adaptable  to  other 
uses  in  addition  to  verifi-cation  of  BACTM  results. 
Since,  regardless  of  the  particular  model  used,  a 
number  of  auxiliary  variables  were  to  be  calculated 
during  the  simulation,  the  additional  calculation  of 
alpha  and  beta  presented  no  increase  in  effort  or 
computation  time. 
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See  appendix  E  for  the  particular  auxiliary  variables 
calculated  during  the  simulations  and  their 
definitions.  Below  is  the  eighth  order  set  of 
equations  of  motion  used  in  the  numerical  simulation 
phase  [7]. 


u  ■  -9  sine  ♦  vr  -  wq  ♦  ^  Cx 
v  *  g  coseslne  ♦  wp  -  ur  ♦  ^ 
w  *  g  cosecose  ♦  uq  -  vp  ♦  fh 


pq 


6  «  q  cose  -  r  sin# 
e  *  p  ♦  q  tanesin*  ♦  r  tanecos* 


(D.2.14) 

(D.2.15) 

(D.2.16) 


(D.2.17) 
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(D. 2.19) 

(D.2.20) 

(0.2.21) 
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APPENDIX  E 
AUXILIARY  VARIBLES 

In  order  to  fully  utilize  the  computational  effort 
expended  during  the  numerical  simulations,  a  number  of 
auxiliary  variables  were  calculated  in  addition  to  the 
trajectories  of  the  state  variables.  Through  these 
variables  a  better  understanding  of  the  behavior  of  the 
aircraft  can  be  obtained  as  well  as  the  chance  to 
monitor  and  verify  some  possible  departure  or  jump 
prediction  parameters.  Table  E. 1  shows  the  variables 
calculated  as  auxiliary  variables  and  the  formulae  used 
to  calculate  the  variables.  A  discussion  of  the 
significance  of  each  of  the  variables  follows  table 
E.  1. 


TABLE  E.  1 


a  *  tan"*(w/u) 

lE.l) 

B  «  sin‘l(v/V) 

(E.2) 

V  *  /iP  ♦  ♦  w^ 

(E.3) 

♦  •  o  *1n#sece  ♦  r  co$#s«ce 

(E.4) 

rtf 

TURNS  •  /  i  dt 

(E.5) 

0 


f  «  f#  ♦  TURNS  (2*) 


CE.6) 
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wi  »  u  sine  -  (v  sine  ♦  w  cosejcose, 
lf 

WI  dt 

0 

E  *  %  mV?  ♦  >i  u*  I*w 

« 

E  -  mW  ♦  Ixpp  ♦  lyqq  ♦  lzrr  +  (rp  +.pr)lxx 

where 

V  *  (uu  ♦  vv  ♦  w w)/V 
If  the  body-axis  velocity  set  Is  used 


u  *  /p2  ♦  q2  ♦  rz 


• 

~ 

°SPIN  * 

n 

I 

'x.-x  " 
i  Sj 

2 

L"1 

• 

Es  *  *  Ive2/(qSb) 


where 

Iy  *  sin2elx  ♦  (sin2ely  ♦  cos2#I2)cos2e  ♦’ 
Ixzcosesin2e 


i 


(E.7) 

CE.8) 

(E.9) 

CE.IO) 

(E. II) 

(E.12) 

(E.13) 

CE.I4) 

CE.I5) 


when  Ixy=Iyz=0 

The  wind  axis  variables  (alpha,  beta  and  total 
velocity)  are  calculated  to  provide  direct  verification 
of  the  results  of  BACTM  computer  runs.  The  heading 
variable,  ,  is  calculated  so  that  an  aircraft's 
horizontal  trajectory  during  a  spin  can  be  calculated. 
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It  is 

also 

used 

in 

the  calculation  of  the 

variable. 

Turns, 

which 

is 

the 

number  of  rotations 

that  have 

occurred  during  a  spin.  This  value  is  of  interest 
since  spin  recovery  techniques  are  often  evaluated  by 
the  number  of  turns  that  occur  between  the  setting  of 
the  spin  recovery  controls  and  the  actual  recovery  from 
the  spin.  Aircraft  vertical  velocity,  WI,  and 
altitude,  h,  are  calculated  to  allow  altitude  to  vary 
during  the  simulation  if  desired  and  to  determine  sink 
rates  during  stalled  and  near  stall  flight.  Vehicle 
kinetic  energy,  E,  and  rate  of  change  of  kinetic 
energy,  E,  are  calculated  to  help  develop  possible 
techniques  to  monitor  and  predict  nonlinear  jump 
behavior.  Dspin  is  a  measurement  of  the  distance  that 
the  current  state  is  from  a  known  spin  condition,  and 
spin  energy,  Es,  is  a  measurement  of  the  kinetic  energy 
due  to  the  spin  motion  of  the  aircraft  alone.  These 
last  two  terms  are  used  to  measure  the  relative 
effectiveness  of  different  spin  recovery  control 
strategies. 
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Fig.  7. 2e :  Side  Slip  Angle  vs.  Aileron 
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Fig.  7.5e*  Side  Slip  Angle  vs.  Rudder 
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